home *** CD-ROM | disk | FTP | other *** search
/ AmigActive 26 / AACD 26.iso / AACD / Graphics / sKulpt / skulpt-src / win-src / Triangulator.cpp < prev    next >
Encoding:
Text File  |  2000-06-12  |  274.5 KB  |  5,898 lines
  1. /*****************************************************************************/
  2. /*                                                                           */
  3. /*      888888888        ,o,                          / 888                  */
  4. /*         888    88o88o  "    o8888o  88o8888o o88888o 888  o88888o         */
  5. /*         888    888    888       88b 888  888 888 888 888 d888  88b        */
  6. /*         888    888    888  o88^o888 888  888 "88888" 888 8888oo888        */
  7. /*         888    888    888 C888  888 888  888  /      888 q888             */
  8. /*         888    888    888  "88o^888 888  888 Cb      888  "88oooo"        */
  9. /*                                              "8oo8D                       */
  10. /*                                                                           */
  11. /*  A Two-Dimensional Quality Mesh Generator and Delaunay Triangulator.      */
  12. /*  (triangle.c)                                                             */
  13. /*                                                                           */
  14. /*  Version 1.3                                                              */
  15. /*  July 19, 1996                                                            */
  16. /*                                                                           */
  17. /*  Copyright 1996                                                           */
  18. /*  Jonathan Richard Shewchuk                                                */
  19. /*  School of Computer Science                                               */
  20. /*  Carnegie Mellon University                                               */
  21. /*  5000 Forbes Avenue                                                       */
  22. /*  Pittsburgh, Pennsylvania  15213-3891                                     */
  23. /*  jrs@cs.cmu.edu                                                           */
  24. /* Severely stripped by SG to reduce size ************************************/
  25. /*                                                                           */
  26. /*  This program may be freely redistributed under the condition that the    */
  27. /*    copyright notices (including this entire header and the copyright      */
  28. /*    notice printed when the `-h' switch is selected) are not removed, and  */
  29. /*    no compensation is received.  Private, research, and institutional     */
  30. /*    use is free.  You may distribute modified versions of this code UNDER  */
  31. /*    THE CONDITION THAT THIS CODE AND ANY MODIFICATIONS MADE TO IT IN THE   */
  32. /*    SAME FILE REMAIN UNDER COPYRIGHT OF THE ORIGINAL AUTHOR, BOTH SOURCE   */
  33. /*    AND OBJECT CODE ARE MADE FREELY AVAILABLE WITHOUT CHARGE, AND CLEAR    */
  34. /*    NOTICE IS GIVEN OF THE MODIFICATIONS.  Distribution of this code as    */
  35. /*    part of a commercial system is permissible ONLY BY DIRECT ARRANGEMENT  */
  36. /*    WITH THE AUTHOR.  (If you are not directly supplying this code to a    */
  37. /*    customer, and you are instead telling them how they can obtain it for  */
  38. /*    free, then you are not required to make any arrangement with me.)      */
  39. /*                                                                           */
  40. /*  Hypertext instructions for Triangle are available on the Web at          */
  41. /*                                                                           */
  42. /*      http://www.cs.cmu.edu/~quake/triangle.html                           */
  43. /*                                                                           */
  44. /*  Some of the references listed below are marked [*].  These are available */
  45. /*    for downloading from the Web page                                      */
  46. /*                                                                           */
  47. /*      http://www.cs.cmu.edu/~quake/triangle.research.html                  */
  48. /*                                                                           */
  49. /*  A paper discussing some aspects of Triangle is available.  See Jonathan  */
  50. /*    Richard Shewchuk, "Triangle:  Engineering a 2D Quality Mesh Generator  */
  51. /*    and Delaunay Triangulator," First Workshop on Applied Computational    */
  52. /*    Geometry, ACM, May 1996.  [*]                                          */
  53. /*                                                                           */
  54. /*  Triangle was created as part of the Archimedes project in the School of  */
  55. /*    Computer Science at Carnegie Mellon University.  Archimedes is a       */
  56. /*    system for compiling parallel finite element solvers.  For further     */
  57. /*    information, see Anja Feldmann, Omar Ghattas, John R. Gilbert, Gary L. */
  58. /*    Miller, David R. O'Hallaron, Eric J. Schwabe, Jonathan R. Shewchuk,    */
  59. /*    and Shang-Hua Teng, "Automated Parallel Solution of Unstructured PDE   */
  60. /*    Problems."  To appear in Communications of the ACM, we hope.           */
  61. /*                                                                           */
  62. /*  The quality mesh generation algorithm is due to Jim Ruppert, "A          */
  63. /*    Delaunay Refinement Algorithm for Quality 2-Dimensional Mesh           */
  64. /*    Generation," Journal of Algorithms 18(3):548-585, May 1995.  [*]       */
  65. /*                                                                           */
  66. /*  My implementation of the divide-and-conquer and incremental Delaunay     */
  67. /*    triangulation algorithms follows closely the presentation of Guibas    */
  68. /*    and Stolfi, even though I use a triangle-based data structure instead  */
  69. /*    of their quad-edge data structure.  (In fact, I originally implemented */
  70. /*    Triangle using the quad-edge data structure, but switching to a        */
  71. /*    triangle-based data structure sped Triangle by a factor of two.)  The  */
  72. /*    mesh manipulation primitives and the two aforementioned Delaunay       */
  73. /*    triangulation algorithms are described by Leonidas J. Guibas and Jorge */
  74. /*    Stolfi, "Primitives for the Manipulation of General Subdivisions and   */
  75. /*    the Computation of Voronoi Diagrams," ACM Transactions on Graphics     */
  76. /*    4(2):74-123, April 1985.                                               */
  77. /*                                                                           */
  78. /*  Their O(n log n) divide-and-conquer algorithm is adapted from Der-Tsai   */
  79. /*    Lee and Bruce J. Schachter, "Two Algorithms for Constructing the       */
  80. /*    Delaunay Triangulation," International Journal of Computer and         */
  81. /*    Information Science 9(3):219-242, 1980.  The idea to improve the       */
  82. /*    divide-and-conquer algorithm by alternating between vertical and       */
  83. /*    horizontal cuts was introduced by Rex A. Dwyer, "A Faster Divide-and-  */
  84. /*    Conquer Algorithm for Constructing Delaunay Triangulations,"           */
  85. /*    Algorithmica 2(2):137-151, 1987.                                       */
  86. /*                                                                           */
  87. /*  The incremental insertion algorithm was first proposed by C. L. Lawson,  */
  88. /*    "Software for C1 Surface Interpolation," in Mathematical Software III, */
  89. /*    John R. Rice, editor, Academic Press, New York, pp. 161-194, 1977.     */
  90. /*    For point location, I use the algorithm of Ernst P. Mucke, Isaac       */
  91. /*    Saias, and Binhai Zhu, "Fast Randomized Point Location Without         */
  92. /*    Preprocessing in Two- and Three-dimensional Delaunay Triangulations,"  */
  93. /*    Proceedings of the Twelfth Annual Symposium on Computational Geometry, */
  94. /*    ACM, May 1996.  [*]  If I were to randomize the order of point         */
  95. /*    insertion (I currently don't bother), their result combined with the   */
  96. /*    result of Leonidas J. Guibas, Donald E. Knuth, and Micha Sharir,       */
  97. /*    "Randomized Incremental Construction of Delaunay and Voronoi           */
  98. /*    Diagrams," Algorithmica 7(4):381-413, 1992, would yield an expected    */
  99. /*    O(n^{4/3}) bound on running time.                                      */
  100. /*                                                                           */
  101. /*  The O(n log n) sweepline Delaunay triangulation algorithm is taken from  */
  102. /*    Steven Fortune, "A Sweepline Algorithm for Voronoi Diagrams",          */
  103. /*    Algorithmica 2(2):153-174, 1987.  A random sample of edges on the      */
  104. /*    boundary of the triangulation are maintained in a splay tree for the   */
  105. /*    purpose of point location.  Splay trees are described by Daniel        */
  106. /*    Dominic Sleator and Robert Endre Tarjan, "Self-Adjusting Binary Search */
  107. /*    Trees," Journal of the ACM 32(3):652-686, July 1985.                   */
  108. /*                                                                           */
  109. /*  The algorithms for exact computation of the signs of determinants are    */
  110. /*    described in Jonathan Richard Shewchuk, "Adaptive Precision Floating-  */
  111. /*    Point Arithmetic and Fast Robust Geometric Predicates," Technical      */
  112. /*    Report CMU-CS-96-140, School of Computer Science, Carnegie Mellon      */
  113. /*    University, Pittsburgh, Pennsylvania, May 1996.  [*]  (Submitted to    */
  114. /*    Discrete & Computational Geometry.)  An abbreviated version appears as */
  115. /*    Jonathan Richard Shewchuk, "Robust Adaptive Floating-Point Geometric   */
  116. /*    Predicates," Proceedings of the Twelfth Annual Symposium on Computa-   */
  117. /*    tional Geometry, ACM, May 1996.  [*]  Many of the ideas for my exact   */
  118. /*    arithmetic routines originate with Douglas M. Priest, "Algorithms for  */
  119. /*    Arbitrary Precision Floating Point Arithmetic," Tenth Symposium on     */
  120. /*    Computer Arithmetic, 132-143, IEEE Computer Society Press, 1991.  [*]  */
  121. /*    Many of the ideas for the correct evaluation of the signs of           */
  122. /*    determinants are taken from Steven Fortune and Christopher J. Van Wyk, */
  123. /*    "Efficient Exact Arithmetic for Computational Geometry," Proceedings   */
  124. /*    of the Ninth Annual Symposium on Computational Geometry, ACM,          */
  125. /*    pp. 163-172, May 1993, and from Steven Fortune, "Numerical Stability   */
  126. /*    of Algorithms for 2D Delaunay Triangulations," International Journal   */
  127. /*    of Computational Geometry & Applications 5(1-2):193-213, March-June    */
  128. /*    1995.                                                                  */
  129. /*                                                                           */
  130. /*  For definitions of and results involving Delaunay triangulations,        */
  131. /*    constrained and conforming versions thereof, and other aspects of      */
  132. /*    triangular mesh generation, see the excellent survey by Marshall Bern  */
  133. /*    and David Eppstein, "Mesh Generation and Optimal Triangulation," in    */
  134. /*    Computing and Euclidean Geometry, Ding-Zhu Du and Frank Hwang,         */
  135. /*    editors, World Scientific, Singapore, pp. 23-90, 1992.                 */
  136. /*                                                                           */
  137. /*  The time for incrementally adding PSLG (planar straight line graph)      */
  138. /*    segments to create a constrained Delaunay triangulation is probably    */
  139. /*    O(n^2) per segment in the worst case and O(n) per edge in the common   */
  140. /*    case, where n is the number of triangles that intersect the segment    */
  141. /*    before it is inserted.  This doesn't count point location, which can   */
  142. /*    be much more expensive.  (This note does not apply to conforming       */
  143. /*    Delaunay triangulations, for which a different method is used to       */
  144. /*    insert segments.)                                                      */
  145. /*                                                                           */
  146. /*  The time for adding segments to a conforming Delaunay triangulation is   */
  147. /*    not clear, but does not depend upon n alone.  In some cases, very      */
  148. /*    small features (like a point lying next to a segment) can cause a      */
  149. /*    single segment to be split an arbitrary number of times.  Of course,   */
  150. /*    floating-point precision is a practical barrier to how much this can   */
  151. /*    happen.                                                                */
  152. /*                                                                           */
  153. /*  The time for deleting a point from a Delaunay triangulation is O(n^2) in */
  154. /*    the worst case and O(n) in the common case, where n is the degree of   */
  155. /*    the point being deleted.  I could improve this to expected O(n) time   */
  156. /*    by "inserting" the neighboring vertices in random order, but n is      */
  157. /*    usually quite small, so it's not worth the bother.  (The O(n) time     */
  158. /*    for random insertion follows from L. Paul Chew, "Building Voronoi      */
  159. /*    Diagrams for Convex Polygons in Linear Expected Time," Technical       */
  160. /*    Report PCS-TR90-147, Department of Mathematics and Computer Science,   */
  161. /*    Dartmouth College, 1990.                                               */
  162. /*                                                                           */
  163. /*  Ruppert's Delaunay refinement algorithm typically generates triangles    */
  164. /*    at a linear rate (constant time per triangle) after the initial        */
  165. /*    triangulation is formed.  There may be pathological cases where more   */
  166. /*    time is required, but these never arise in practice.                   */
  167. /*                                                                           */
  168. /*  The segment intersection formulae are straightforward.  If you want to   */
  169. /*    see them derived, see Franklin Antonio.  "Faster Line Segment          */
  170. /*    Intersection."  In Graphics Gems III (David Kirk, editor), pp. 199-    */
  171. /*    202.  Academic Press, Boston, 1992.                                    */
  172. /*                                                                           */
  173. /*  If you make any improvements to this code, please please please let me   */
  174. /*    know, so that I may obtain the improvements.  Even if you don't change */
  175. /*    the code, I'd still love to hear what it's being used for.             */
  176. /*                                                                           */
  177. /*  Disclaimer:  Neither I nor Carnegie Mellon warrant this code in any way  */
  178. /*    whatsoever.  This code is provided "as-is".  Use at your own risk.     */
  179. /*                                                                           */
  180. /*****************************************************************************/
  181.  
  182. /* On some machines, the exact arithmetic routines might be defeated by the  */
  183. /*   use of internal extended precision floating-point registers.  Sometimes */
  184. /*   this problem can be fixed by defining certain values to be volatile,    */
  185. /*   thus forcing them to be stored to memory and rounded off.  This isn't   */
  186. /*   a great solution, though, as it slows Triangle down.                    */
  187. /*                                                                           */
  188. /* To try this out, write "#define INEXACT volatile" below.  Normally,       */
  189. /*   however, INEXACT should be defined to be nothing.  ("#define INEXACT".) */
  190.  
  191. #define INEXACT
  192.  
  193. /* For efficiency, a variety of data structures are allocated in bulk.  The  */
  194. /*   following constants determine how many of each structure is allocated   */
  195. /*   at once.                                                                */
  196.  
  197. #define TRIPERBLOCK 4092           /* Number of triangles allocated at once. */
  198. #define SHELLEPERBLOCK 508       /* Number of shell edges allocated at once. */
  199. #define POINTPERBLOCK 4092            /* Number of points allocated at once. */
  200. #define VIRUSPERBLOCK 1020   /* Number of virus triangles allocated at once. */
  201.  
  202. /* The point marker DEADPOINT is an arbitrary number chosen large enough to  */
  203. /*   (hopefully) not conflict with user boundary markers.  Make sure that it */
  204. /*   is small enough to fit into your machine's integer size.                */
  205. #define DEADPOINT -1073741824
  206.  
  207. /* Used for the point location scheme of Mucke, Saias, and Zhu, to decide    */
  208. /*   how large a random sample of triangles to inspect.                      */
  209. #define SAMPLEFACTOR 11
  210.  
  211. /* A number that speaks for itself, every kissable digit.                    */
  212. #define PI 3.141592653589793238462643383279502884197169399375105820974944592308
  213.  
  214. // Patch SG pour intégration MSVC6 / projet sKulpt
  215. #define STRICT
  216. extern void    vTrace(char *Str, ...);
  217. #include <stdlib.h>
  218. // End patch
  219.  
  220. #include <stdio.h>
  221. #include <string.h>
  222. #include <math.h>
  223. #include "triangulator.h"
  224.  
  225. /* A few forward declarations.                                               */
  226. void poolrestart(struct memorypool *pool);
  227.  
  228. /* Labels that signify whether a record consists primarily of pointers or of */
  229. /*   floating-point words.  Used to make decisions about data alignment.     */
  230.  
  231. enum wordtype {POINTER, FLOATINGPOINT};
  232.  
  233. /* Labels that signify the result of point location.  The result of a        */
  234. /*   search indicates that the point falls in the interior of a triangle, on */
  235. /*   an edge, on a vertex, or outside the mesh.                              */
  236.  
  237. enum locateresult {INTRIANGLE, ONEDGE, ONVERTEX, OUTSIDE};
  238.  
  239. /* Labels that signify the result of site insertion.  The result indicates   */
  240. /*   that the point was inserted with complete success, was inserted but     */
  241. /*   encroaches on a segment, was not inserted because it lies on a segment, */
  242. /*   or was not inserted because another point occupies the same location.   */
  243.  
  244. enum insertsiteresult {SUCCESSFULPOINT, VIOLATINGPOINT, DUPLICATEPOINT};
  245.  
  246. /* Labels that signify the result of direction finding.  The result          */
  247. /*   indicates that a segment connecting the two query points falls within   */
  248. /*   the direction triangle, along the left edge of the direction triangle,  */
  249. /*   or along the right edge of the direction triangle.                      */
  250.  
  251. enum finddirectionresult {WITHIN, LEFTCOLLINEAR, RIGHTCOLLINEAR};
  252.  
  253. /*****************************************************************************/
  254. /*                                                                           */
  255. /*  The basic mesh data structures                                           */
  256. /*                                                                           */
  257. /*  There are three:  points, triangles, and shell edges (abbreviated        */
  258. /*  `shelle').  These three data structures, linked by pointers, comprise    */
  259. /*  the mesh.  A point simply represents a point in space and its properties.*/
  260. /*  A triangle is a triangle.  A shell edge is a special data structure used */
  261. /*  to represent impenetrable segments in the mesh (including the outer      */
  262. /*  boundary, boundaries of holes, and internal boundaries separating two    */
  263. /*  triangulated regions).  Shell edges represent boundaries defined by the  */
  264. /*  user that triangles may not lie across.                                  */
  265. /*                                                                           */
  266. /*  A triangle consists of a list of three vertices, a list of three         */
  267. /*  adjoining triangles, a list of three adjoining shell edges (when shell   */
  268. /*  edges are used), an arbitrary number of optional user-defined floating-  */
  269. /*  point attributes, and an optional area constraint.  The latter is an     */
  270. /*  upper bound on the permissible area of each triangle in a region, used   */
  271. /*  for mesh refinement.                                                     */
  272. /*                                                                           */
  273. /*  For a triangle on a boundary of the mesh, some or all of the neighboring */
  274. /*  triangles may not be present.  For a triangle in the interior of the     */
  275. /*  mesh, often no neighboring shell edges are present.  Such absent         */
  276. /*  triangles and shell edges are never represented by NULL pointers; they   */
  277. /*  are represented by two special records:  `dummytri', the triangle that   */
  278. /*  fills "outer space", and `dummysh', the omnipresent shell edge.          */
  279. /*  `dummytri' and `dummysh' are used for several reasons; for instance,     */
  280. /*  they can be dereferenced and their contents examined without causing the */
  281. /*  memory protection exception that would occur if NULL were dereferenced.  */
  282. /*                                                                           */
  283. /*  However, it is important to understand that a triangle includes other    */
  284. /*  information as well.  The pointers to adjoining vertices, triangles, and */
  285. /*  shell edges are ordered in a way that indicates their geometric relation */
  286. /*  to each other.  Furthermore, each of these pointers contains orientation */
  287. /*  information.  Each pointer to an adjoining triangle indicates which face */
  288. /*  of that triangle is contacted.  Similarly, each pointer to an adjoining  */
  289. /*  shell edge indicates which side of that shell edge is contacted, and how */
  290. /*  the shell edge is oriented relative to the triangle.                     */
  291. /*                                                                           */
  292. /*  Shell edges are found abutting edges of triangles; either sandwiched     */
  293. /*  between two triangles, or resting against one triangle on an exterior    */
  294. /*  boundary or hole boundary.                                               */
  295. /*                                                                           */
  296. /*  A shell edge consists of a list of two vertices, a list of two           */
  297. /*  adjoining shell edges, and a list of two adjoining triangles.  One of    */
  298. /*  the two adjoining triangles may not be present (though there should      */
  299. /*  always be one), and neighboring shell edges might not be present.        */
  300. /*  Shell edges also store a user-defined integer "boundary marker".         */
  301. /*  Typically, this integer is used to indicate what sort of boundary        */
  302. /*  conditions are to be applied at that location in a finite element        */
  303. /*  simulation.                                                              */
  304. /*                                                                           */
  305. /*  Like triangles, shell edges maintain information about the relative      */
  306. /*  orientation of neighboring objects.                                      */
  307. /*                                                                           */
  308. /*  Points are relatively simple.  A point is a list of floating point       */
  309. /*  numbers, starting with the x, and y coordinates, followed by an          */
  310. /*  arbitrary number of optional user-defined floating-point attributes,     */
  311. /*  followed by an integer boundary marker.  During the segment insertion    */
  312. /*  phase, there is also a pointer from each point to a triangle that may    */
  313. /*  contain it.  Each pointer is not always correct, but when one is, it     */
  314. /*  speeds up segment insertion.  These pointers are assigned values once    */
  315. /*  at the beginning of the segment insertion phase, and are not used or     */
  316. /*  updated at any other time.  Edge swapping during segment insertion will  */
  317. /*  render some of them incorrect.  Hence, don't rely upon them for          */
  318. /*  anything.  For the most part, points do not have any information about   */
  319. /*  what triangles or shell edges they are linked to.                        */
  320. /*                                                                           */
  321. /*****************************************************************************/
  322.  
  323. /*****************************************************************************/
  324. /*                                                                           */
  325. /*  Handles                                                                  */
  326. /*                                                                           */
  327. /*  The oriented triangle (`triedge') and oriented shell edge (`edge') data  */
  328. /*  structures defined below do not themselves store any part of the mesh.   */
  329. /*  The mesh itself is made of `triangle's, `shelle's, and `point's.         */
  330. /*                                                                           */
  331. /*  Oriented triangles and oriented shell edges will usually be referred to  */
  332. /*  as "handles".  A handle is essentially a pointer into the mesh; it       */
  333. /*  allows you to "hold" one particular part of the mesh.  Handles are used  */
  334. /*  to specify the regions in which one is traversing and modifying the mesh.*/
  335. /*  A single `triangle' may be held by many handles, or none at all.  (The   */
  336. /*  latter case is not a memory leak, because the triangle is still          */
  337. /*  connected to other triangles in the mesh.)                               */
  338. /*                                                                           */
  339. /*  A `triedge' is a handle that holds a triangle.  It holds a specific side */
  340. /*  of the triangle.  An `edge' is a handle that holds a shell edge.  It     */
  341. /*  holds either the left or right side of the edge.                         */
  342. /*                                                                           */
  343. /*  Navigation about the mesh is accomplished through a set of mesh          */
  344. /*  manipulation primitives, further below.  Many of these primitives take   */
  345. /*  a handle and produce a new handle that holds the mesh near the first     */
  346. /*  handle.  Other primitives take two handles and glue the corresponding    */
  347. /*  parts of the mesh together.  The exact position of the handles is        */
  348. /*  important.  For instance, when two triangles are glued together by the   */
  349. /*  bond() primitive, they are glued by the sides on which the handles lie.  */
  350. /*                                                                           */
  351. /*  Because points have no information about which triangles they are        */
  352. /*  attached to, I commonly represent a point by use of a handle whose       */
  353. /*  origin is the point.  A single handle can simultaneously represent a     */
  354. /*  triangle, an edge, and a point.                                          */
  355. /*                                                                           */
  356. /*****************************************************************************/
  357.  
  358. /* The triangle data structure.  Each triangle contains three pointers to    */
  359. /*   adjoining triangles, plus three pointers to vertex points, plus three   */
  360. /*   pointers to shell edges (defined below; these pointers are usually      */
  361. /*   `dummysh').  It may or may not also contain user-defined attributes     */
  362. /*   and/or a floating-point "area constraint".  It may also contain extra   */
  363. /*   pointers for nodes, when the user asks for high-order elements.         */
  364. /*   Because the size and structure of a `triangle' is not decided until     */
  365. /*   runtime, I haven't simply defined the type `triangle' to be a struct.   */
  366.  
  367. typedef double **triangle;            /* Really:  typedef triangle *triangle   */
  368.  
  369. /* An oriented triangle:  includes a pointer to a triangle and orientation.  */
  370. /*   The orientation denotes an edge of the triangle.  Hence, there are      */
  371. /*   three possible orientations.  By convention, each edge is always        */
  372. /*   directed to point counterclockwise about the corresponding triangle.    */
  373.  
  374. struct triedge {
  375.   triangle *tri;
  376.   int orient;                                         /* Ranges from 0 to 2. */
  377. };
  378.  
  379. /* The shell data structure.  Each shell edge contains two pointers to       */
  380. /*   adjoining shell edges, plus two pointers to vertex points, plus two     */
  381. /*   pointers to adjoining triangles, plus one shell marker.                 */
  382.  
  383. typedef double **shelle;                  /* Really:  typedef shelle *shelle   */
  384.  
  385. /* An oriented shell edge:  includes a pointer to a shell edge and an        */
  386. /*   orientation.  The orientation denotes a side of the edge.  Hence, there */
  387. /*   are two possible orientations.  By convention, the edge is always       */
  388. /*   directed so that the "side" denoted is the right side of the edge.      */
  389.  
  390. struct edge {
  391.   shelle *sh;
  392.   int shorient;                                       /* Ranges from 0 to 1. */
  393. };
  394.  
  395. /* The point data structure.  Each point is actually an array of REALs.      */
  396. /*   The number of REALs is unknown until runtime.  An integer boundary      */
  397. /*   marker, and sometimes a pointer to a triangle, is appended after the    */
  398. /*   REALs.                                                                  */
  399.  
  400. typedef double *point;
  401.  
  402. /* A type used to allocate memory.  firstblock is the first block of items.  */
  403. /*   nowblock is the block from which items are currently being allocated.   */
  404. /*   nextitem points to the next slab of free memory for an item.            */
  405. /*   deaditemstack is the head of a linked list (stack) of deallocated items */
  406. /*   that can be recycled.  unallocateditems is the number of items that     */
  407. /*   remain to be allocated from nowblock.                                   */
  408. /*                                                                           */
  409. /* Traversal is the process of walking through the entire list of items, and */
  410. /*   is separate from allocation.  Note that a traversal will visit items on */
  411. /*   the "deaditemstack" stack as well as live items.  pathblock points to   */
  412. /*   the block currently being traversed.  pathitem points to the next item  */
  413. /*   to be traversed.  pathitemsleft is the number of items that remain to   */
  414. /*   be traversed in pathblock.                                              */
  415. /*                                                                           */
  416. /* itemwordtype is set to POINTER or FLOATINGPOINT, and is used to suggest   */
  417. /*   what sort of word the record is primarily made up of.  alignbytes       */
  418. /*   determines how new records should be aligned in memory.  itembytes and  */
  419. /*   itemwords are the length of a record in bytes (after rounding up) and   */
  420. /*   words.  itemsperblock is the number of items allocated at once in a     */
  421. /*   single block.  items is the number of currently allocated items.        */
  422. /*   maxitems is the maximum number of items that have been allocated at     */
  423. /*   once; it is the current number of items plus the number of records kept */
  424. /*   on deaditemstack.                                                       */
  425.  
  426. struct memorypool {
  427.   void **firstblock, **nowblock;
  428.   void *nextitem;
  429.   void *deaditemstack;
  430.   void **pathblock;
  431.   void *pathitem;
  432.   enum wordtype itemwordtype;
  433.   int alignbytes;
  434.   int itembytes, itemwords;
  435.   int itemsperblock;
  436.   long items, maxitems;
  437.   int unallocateditems;
  438.   int pathitemsleft;
  439. };
  440.  
  441. /* Variables used to allocate memory for triangles, shell edges, points,     */
  442. /*   viri (triangles being eaten), bad (encroached) segments, bad (skinny    */
  443. /*   or too large) triangles, and splay tree nodes.                          */
  444.  
  445. struct memorypool triangles;
  446. struct memorypool shelles;
  447. struct memorypool points;
  448. struct memorypool viri;
  449. struct memorypool badsegments;
  450. struct memorypool badtriangles;
  451. struct memorypool splaynodes;
  452.  
  453. /* Variables that maintain the bad triangle queues.  The tails are pointers  */
  454. /*   to the pointers that have to be filled in to enqueue an item.           */
  455.  
  456. double xmin, xmax, ymin, ymax;                              /* x and y bounds. */
  457. double xminextreme;        /* Nonexistent x value used as a flag in sweepline. */
  458. int inpoints;                                     /* Number of input points. */
  459. int insegments;                                 /* Number of input segments. */
  460. int holes;                                         /* Number of input holes. */
  461. int regions;                                     /* Number of input regions. */
  462. long edges;                                       /* Number of output edges. */
  463. int mesh_dim;                                  /* Dimension (ought to be 2). */
  464. int nextras;                              /* Number of attributes per point. */
  465. int eextras;                           /* Number of attributes per triangle. */
  466. long hullsize;                            /* Number of edges of convex hull. */
  467. int triwords;                                   /* Total words per triangle. */
  468. int shwords;                                  /* Total words per shell edge. */
  469. int pointmarkindex;             /* Index to find boundary marker of a point. */
  470. int point2triindex;         /* Index to find a triangle adjacent to a point. */
  471. int highorderindex;    /* Index to find extra nodes for high-order elements. */
  472. int elemattribindex;              /* Index to find attributes of a triangle. */
  473. int areaboundindex;               /* Index to find area bound of a triangle. */
  474. int checksegments;           /* Are there segments in the triangulation yet? */
  475. long samples;                /* Number of random samples for point location. */
  476. unsigned long randomseed;                     /* Current random number seed. */
  477.  
  478. double splitter;       /* Used to split double factors for exact multiplication. */
  479. double epsilon;                             /* Floating-point machine epsilon. */
  480. double resulterrbound;
  481. double ccwerrboundA, ccwerrboundB, ccwerrboundC;
  482. double iccerrboundA, iccerrboundB, iccerrboundC;
  483.  
  484. long incirclecount;                   /* Number of incircle tests performed. */
  485. long counterclockcount;       /* Number of counterclockwise tests performed. */
  486. long hyperbolacount;        /* Number of right-of-hyperbola tests performed. */
  487. long circumcentercount;    /* Number of circumcenter calculations performed. */
  488. long circletopcount;         /* Number of circle top calculations performed. */
  489.  
  490. /* Switches for the triangulator.                                            */
  491. /*   useshelles: -p, -r, -q, or -c switch; determines whether shell edges    */
  492. /*     are used at all.                                                      */
  493. int useshelles;
  494. int order;
  495. double minangle, goodangle;
  496. double maxarea;
  497.  
  498. /* Triangular bounding box points.                                           */
  499.  
  500. point infpoint1, infpoint2, infpoint3;
  501.  
  502. /* Pointer to the `triangle' that occupies all of "outer space".             */
  503.  
  504. triangle *dummytri;
  505. triangle *dummytribase;      /* Keep base address so we can free() it later. */
  506.  
  507. /* Pointer to the omnipresent shell edge.  Referenced by any triangle or     */
  508. /*   shell edge that isn't really connected to a shell edge at that          */
  509. /*   location.                                                               */
  510.  
  511. shelle *dummysh;
  512. shelle *dummyshbase;         /* Keep base address so we can free() it later. */
  513.  
  514. /* Pointer to a recently visited triangle.  Improves point location if       */
  515. /*   proximate points are inserted sequentially.                             */
  516.  
  517. struct triedge recenttri;
  518.  
  519. /*****************************************************************************/
  520. /*                                                                           */
  521. /*  Mesh manipulation primitives.  Each triangle contains three pointers to  */
  522. /*  other triangles, with orientations.  Each pointer points not to the      */
  523. /*  first byte of a triangle, but to one of the first three bytes of a       */
  524. /*  triangle.  It is necessary to extract both the triangle itself and the   */
  525. /*  orientation.  To save memory, I keep both pieces of information in one   */
  526. /*  pointer.  To make this possible, I assume that all triangles are aligned */
  527. /*  to four-byte boundaries.  The `decode' routine below decodes a pointer,  */
  528. /*  extracting an orientation (in the range 0 to 2) and a pointer to the     */
  529. /*  beginning of a triangle.  The `encode' routine compresses a pointer to a */
  530. /*  triangle and an orientation into a single pointer.  My assumptions that  */
  531. /*  triangles are four-byte-aligned and that the `unsigned long' type is     */
  532. /*  long enough to hold a pointer are two of the few kludges in this program.*/
  533. /*                                                                           */
  534. /*  Shell edges are manipulated similarly.  A pointer to a shell edge        */
  535. /*  carries both an address and an orientation in the range 0 to 1.          */
  536. /*                                                                           */
  537. /*  The other primitives take an oriented triangle or oriented shell edge,   */
  538. /*  and return an oriented triangle or oriented shell edge or point; or they */
  539. /*  change the connections in the data structure.                            */
  540. /*                                                                           */
  541. /*****************************************************************************/
  542.  
  543. /********* Mesh manipulation primitives begin here                   *********/
  544. /**                                                                         **/
  545. /**                                                                         **/
  546.  
  547. /* Fast lookup arrays to speed some of the mesh manipulation primitives.     */
  548.  
  549. int plus1mod3[3] = {1, 2, 0};
  550. int minus1mod3[3] = {2, 0, 1};
  551.  
  552. /********* Primitives for triangles                                  *********/
  553. /*                                                                           */
  554. /*                                                                           */
  555.  
  556. /* decode() converts a pointer to an oriented triangle.  The orientation is  */
  557. /*   extracted from the two least significant bits of the pointer.           */
  558.  
  559. #define decode(ptr, triedge)                                                  \
  560.   (triedge).orient = (int) ((unsigned long) (ptr) & (unsigned long) 3l);      \
  561.   (triedge).tri = (triangle *)                                                \
  562.                   ((unsigned long) (ptr) ^ (unsigned long) (triedge).orient)
  563.  
  564. /* encode() compresses an oriented triangle into a single pointer.  It       */
  565. /*   relies on the assumption that all triangles are aligned to four-byte    */
  566. /*   boundaries, so the two least significant bits of (triedge).tri are zero.*/
  567.  
  568. #define encode(triedge)                                                       \
  569.   (triangle) ((unsigned long) (triedge).tri | (unsigned long) (triedge).orient)
  570.  
  571. /* The following edge manipulation primitives are all described by Guibas    */
  572. /*   and Stolfi.  However, they use an edge-based data structure, whereas I  */
  573. /*   am using a triangle-based data structure.                               */
  574.  
  575. /* sym() finds the abutting triangle, on the same edge.  Note that the       */
  576. /*   edge direction is necessarily reversed, because triangle/edge handles   */
  577. /*   are always directed counterclockwise around the triangle.               */
  578.  
  579. #define sym(triedge1, triedge2)                                               \
  580.   ptr = (triedge1).tri[(triedge1).orient];                                    \
  581.   decode(ptr, triedge2);
  582.  
  583. #define symself(triedge)                                                      \
  584.   ptr = (triedge).tri[(triedge).orient];                                      \
  585.   decode(ptr, triedge);
  586.  
  587. /* lnext() finds the next edge (counterclockwise) of a triangle.             */
  588.  
  589. #define lnext(triedge1, triedge2)                                             \
  590.   (triedge2).tri = (triedge1).tri;                                            \
  591.   (triedge2).orient = plus1mod3[(triedge1).orient]
  592.  
  593. #define lnextself(triedge)                                                    \
  594.   (triedge).orient = plus1mod3[(triedge).orient]
  595.  
  596. /* lprev() finds the previous edge (clockwise) of a triangle.                */
  597.  
  598. #define lprev(triedge1, triedge2)                                             \
  599.   (triedge2).tri = (triedge1).tri;                                            \
  600.   (triedge2).orient = minus1mod3[(triedge1).orient]
  601.  
  602. #define lprevself(triedge)                                                    \
  603.   (triedge).orient = minus1mod3[(triedge).orient]
  604.  
  605. /* onext() spins counterclockwise around a point; that is, it finds the next */
  606. /*   edge with the same origin in the counterclockwise direction.  This edge */
  607. /*   will be part of a different triangle.                                   */
  608.  
  609. #define onext(triedge1, triedge2)                                             \
  610.   lprev(triedge1, triedge2);                                                  \
  611.   symself(triedge2);
  612.  
  613. #define onextself(triedge)                                                    \
  614.   lprevself(triedge);                                                         \
  615.   symself(triedge);
  616.  
  617. /* oprev() spins clockwise around a point; that is, it finds the next edge   */
  618. /*   with the same origin in the clockwise direction.  This edge will be     */
  619. /*   part of a different triangle.                                           */
  620.  
  621. #define oprev(triedge1, triedge2)                                             \
  622.   sym(triedge1, triedge2);                                                    \
  623.   lnextself(triedge2);
  624.  
  625. #define oprevself(triedge)                                                    \
  626.   symself(triedge);                                                           \
  627.   lnextself(triedge);
  628.  
  629. /* These primitives determine or set the origin, destination, or apex of a   */
  630. /* triangle.                                                                 */
  631.  
  632. #define org(triedge, pointptr)                                                \
  633.   pointptr = (point) (triedge).tri[plus1mod3[(triedge).orient] + 3]
  634.  
  635. #define dest(triedge, pointptr)                                               \
  636.   pointptr = (point) (triedge).tri[minus1mod3[(triedge).orient] + 3]
  637.  
  638. #define apex(triedge, pointptr)                                               \
  639.   pointptr = (point) (triedge).tri[(triedge).orient + 3]
  640.  
  641. #define setorg(triedge, pointptr)                                             \
  642.   (triedge).tri[plus1mod3[(triedge).orient] + 3] = (triangle) pointptr
  643.  
  644. #define setdest(triedge, pointptr)                                            \
  645.   (triedge).tri[minus1mod3[(triedge).orient] + 3] = (triangle) pointptr
  646.  
  647. #define setapex(triedge, pointptr)                                            \
  648.   (triedge).tri[(triedge).orient + 3] = (triangle) pointptr
  649.  
  650. /* Bond two triangles together.                                              */
  651.  
  652. #define bond(triedge1, triedge2)                                              \
  653.   (triedge1).tri[(triedge1).orient] = encode(triedge2);                       \
  654.   (triedge2).tri[(triedge2).orient] = encode(triedge1)
  655.  
  656. /* Dissolve a bond (from one side).  Note that the other triangle will still */
  657. /*   think it's connected to this triangle.  Usually, however, the other     */
  658. /*   triangle is being deleted entirely, or bonded to another triangle, so   */
  659. /*   it doesn't matter.                                                      */
  660.  
  661. #define dissolve(triedge)                                                     \
  662.   (triedge).tri[(triedge).orient] = (triangle) dummytri
  663.  
  664. /* Copy a triangle/edge handle.                                              */
  665.  
  666. #define triedgecopy(triedge1, triedge2)                                       \
  667.   (triedge2).tri = (triedge1).tri;                                            \
  668.   (triedge2).orient = (triedge1).orient
  669.  
  670. /* Test for equality of triangle/edge handles.                               */
  671.  
  672. #define triedgeequal(triedge1, triedge2)                                      \
  673.   (((triedge1).tri == (triedge2).tri) &&                                      \
  674.    ((triedge1).orient == (triedge2).orient))
  675.  
  676. /* Primitives to infect or cure a triangle with the virus.  These rely on    */
  677. /*   the assumption that all shell edges are aligned to four-byte boundaries.*/
  678.  
  679. #define infect(triedge)                                                       \
  680.   (triedge).tri[6] = (triangle)                                               \
  681.                      ((unsigned long) (triedge).tri[6] | (unsigned long) 2l)
  682.  
  683. #define uninfect(triedge)                                                     \
  684.   (triedge).tri[6] = (triangle)                                               \
  685.                      ((unsigned long) (triedge).tri[6] & ~ (unsigned long) 2l)
  686.  
  687. /* Test a triangle for viral infection.                                      */
  688.  
  689. #define infected(triedge)                                                     \
  690.   (((unsigned long) (triedge).tri[6] & (unsigned long) 2l) != 0)
  691.  
  692. /* Check or set a triangle's attributes.                                     */
  693.  
  694. #define elemattribute(triedge, attnum)                                        \
  695.   ((double *) (triedge).tri)[elemattribindex + (attnum)]
  696.  
  697. #define setelemattribute(triedge, attnum, value)                              \
  698.   ((double *) (triedge).tri)[elemattribindex + (attnum)] = value
  699.  
  700. /********* Primitives for shell edges                                *********/
  701. /*                                                                           */
  702. /*                                                                           */
  703.  
  704. /* sdecode() converts a pointer to an oriented shell edge.  The orientation  */
  705. /*   is extracted from the least significant bit of the pointer.  The two    */
  706. /*   least significant bits (one for orientation, one for viral infection)   */
  707. /*   are masked out to produce the real pointer.                             */
  708.  
  709. #define sdecode(sptr, edge)                                                   \
  710.   (edge).shorient = (int) ((unsigned long) (sptr) & (unsigned long) 1l);      \
  711.   (edge).sh = (shelle *)                                                      \
  712.               ((unsigned long) (sptr) & ~ (unsigned long) 3l)
  713.  
  714. /* sencode() compresses an oriented shell edge into a single pointer.  It    */
  715. /*   relies on the assumption that all shell edges are aligned to two-byte   */
  716. /*   boundaries, so the least significant bit of (edge).sh is zero.          */
  717.  
  718. #define sencode(edge)                                                         \
  719.   (shelle) ((unsigned long) (edge).sh | (unsigned long) (edge).shorient)
  720.  
  721. /* ssym() toggles the orientation of a shell edge.                           */
  722. #define ssymself(edge)                                                        \
  723.   (edge).shorient = 1 - (edge).shorient
  724.  
  725. /* spivot() finds the other shell edge (from the same segment) that shares   */
  726. /*   the same origin.                                                        */
  727.  
  728. #define spivot(edge1, edge2)                                                  \
  729.   sptr = (edge1).sh[(edge1).shorient];                                        \
  730.   sdecode(sptr, edge2)
  731.  
  732. /* These primitives determine or set the origin or destination of a shell    */
  733. /*   edge.                                                                   */
  734.  
  735. #define setsorg(edge, pointptr)                                               \
  736.   (edge).sh[2 + (edge).shorient] = (shelle) pointptr
  737.  
  738. #define setsdest(edge, pointptr)                                              \
  739.   (edge).sh[3 - (edge).shorient] = (shelle) pointptr
  740.  
  741. /* These primitives read or set a shell marker.  Shell markers are used to   */
  742. /*   hold user boundary information.                                         */
  743.  
  744. #define mark(edge)  (* (int *) ((edge).sh + 6))
  745.  
  746. #define setmark(edge, value)                                                  \
  747.   * (int *) ((edge).sh + 6) = value
  748.  
  749. /* Bond two shell edges together.                                            */
  750.  
  751. #define sbond(edge1, edge2)                                                   \
  752.   (edge1).sh[(edge1).shorient] = sencode(edge2);                              \
  753.   (edge2).sh[(edge2).shorient] = sencode(edge1)
  754.  
  755. /* Copy a shell edge.                                                        */
  756.  
  757. #define shellecopy(edge1, edge2)                                              \
  758.   (edge2).sh = (edge1).sh;                                                    \
  759.   (edge2).shorient = (edge1).shorient
  760.  
  761. /********* Primitives for interacting triangles and shell edges      *********/
  762. /*                                                                           */
  763. /*                                                                           */
  764.  
  765. /* tspivot() finds a shell edge abutting a triangle.                         */
  766.  
  767. #define tspivot(triedge, edge)                                                \
  768.   sptr = (shelle) (triedge).tri[6 + (triedge).orient];                        \
  769.   sdecode(sptr, edge)
  770.  
  771. /* Bond a triangle to a shell edge.                                          */
  772.  
  773. #define tsbond(triedge, edge)                                                 \
  774.   (triedge).tri[6 + (triedge).orient] = (triangle) sencode(edge);             \
  775.   (edge).sh[4 + (edge).shorient] = (shelle) encode(triedge)
  776.  
  777. /* Dissolve a bond (from the triangle side).                                 */
  778.  
  779. #define tsdissolve(triedge)                                                   \
  780.   (triedge).tri[6 + (triedge).orient] = (triangle) dummysh
  781.  
  782. /* Dissolve a bond (from the shell edge side).                               */
  783.  
  784. #define stdissolve(edge)                                                      \
  785.   (edge).sh[4 + (edge).shorient] = (shelle) dummytri
  786.  
  787. /********* Primitives for points                                     *********/
  788. /*                                                                           */
  789. /*                                                                           */
  790.  
  791. #define pointmark(pt)  ((int *) (pt))[pointmarkindex]
  792.  
  793. #define setpointmark(pt, value)                                               \
  794.   ((int *) (pt))[pointmarkindex] = value
  795.  
  796. #define point2tri(pt)  ((triangle *) (pt))[point2triindex]
  797.  
  798. #define setpoint2tri(pt, value)                                               \
  799.   ((triangle *) (pt))[point2triindex] = value
  800.  
  801. /**                                                                         **/
  802. /**                                                                         **/
  803. /********* Mesh manipulation primitives end here                     *********/
  804.  
  805. /********* User interaction routines begin here                      *********/
  806. /**                                                                         **/
  807. /**                                                                         **/
  808.  
  809. /*****************************************************************************/
  810. /*                                                                           */
  811. /*  internalerror()   Ask the user to send me the defective product.  Exit.  */
  812. /*                                                                           */
  813. /*****************************************************************************/
  814.  
  815. void internalerror(void)
  816. {
  817.   vTrace("*** E0031 : Erreur interne. Communiquez le contexte à stephane.guillard@steria.fr");
  818.   exit(1);
  819. }
  820.  
  821. /*****************************************************************************/
  822. /*                                                                           */
  823. /*  parsecommandline()   Read the command line, identify switches, and set   */
  824. /*                       up options and file names.                          */
  825. /*                                                                           */
  826. /*  The effects of this routine are felt entirely through global variables.  */
  827. /*                                                                           */
  828. /*****************************************************************************/
  829.  
  830. void parsecommandline(int argc, char **argv)
  831. {
  832.   order = 1;
  833.   minangle = 0.0;
  834.   maxarea = -1.0;
  835.   
  836.   useshelles = 1;
  837.   goodangle = cos(minangle * PI / 180.0);
  838.   goodangle *= goodangle;
  839. }
  840.  
  841. /**                                                                         **/
  842. /**                                                                         **/
  843. /********* User interaction routines begin here                      *********/
  844.  
  845. /********* Memory management routines begin here                     *********/
  846. /**                                                                         **/
  847. /**                                                                         **/
  848.  
  849. /*****************************************************************************/
  850. /*                                                                           */
  851. /*  poolinit()   Initialize a pool of memory for allocation of items.        */
  852. /*                                                                           */
  853. /*  This routine initializes the machinery for allocating items.  A `pool'   */
  854. /*  is created whose records have size at least `bytecount'.  Items will be  */
  855. /*  allocated in `itemcount'-item blocks.  Each item is assumed to be a      */
  856. /*  collection of words, and either pointers or floating-point values are    */
  857. /*  assumed to be the "primary" word type.  (The "primary" word type is used */
  858. /*  to determine alignment of items.)  If `alignment' isn't zero, all items  */
  859. /*  will be `alignment'-byte aligned in memory.  `alignment' must be either  */
  860. /*  a multiple or a factor of the primary word size; powers of two are safe. */
  861. /*  `alignment' is normally used to create a few unused bits at the bottom   */
  862. /*  of each item's pointer, in which information may be stored.              */
  863. /*                                                                           */
  864. /*  Don't change this routine unless you understand it.                      */
  865. /*                                                                           */
  866. /*****************************************************************************/
  867.  
  868. void poolinit(struct memorypool *pool,int  bytecount, int itemcount, enum wordtype wtype, int alignment)
  869. {
  870.   int wordsize;
  871.  
  872.   /* Initialize values in the pool. */
  873.   pool->itemwordtype = wtype;
  874.   wordsize = (pool->itemwordtype == POINTER) ? sizeof(void *) : sizeof(double);
  875.   /* Find the proper alignment, which must be at least as large as:   */
  876.   /*   - The parameter `alignment'.                                   */
  877.   /*   - The primary word type, to avoid unaligned accesses.          */
  878.   /*   - sizeof(void *), so the stack of dead items can be maintained */
  879.   /*       without unaligned accesses.                                */
  880.   if (alignment > wordsize) {
  881.     pool->alignbytes = alignment;
  882.   } else {
  883.     pool->alignbytes = wordsize;
  884.   }
  885.   if (sizeof(void *) > pool->alignbytes) {
  886.     pool->alignbytes = sizeof(void *);
  887.   }
  888.   pool->itemwords = ((bytecount + pool->alignbytes - 1) / pool->alignbytes)
  889.                   * (pool->alignbytes / wordsize);
  890.   pool->itembytes = pool->itemwords * wordsize;
  891.   pool->itemsperblock = itemcount;
  892.  
  893.   /* Allocate a block of items.  Space for `itemsperblock' items and one    */
  894.   /*   pointer (to point to the next block) are allocated, as well as space */
  895.   /*   to ensure alignment of the items.                                    */
  896.   pool->firstblock = (void **) malloc(pool->itemsperblock * pool->itembytes
  897.                                       + sizeof(void *) + pool->alignbytes);
  898.   if (pool->firstblock == (void **) NULL) {
  899.     vTrace("Error:  Out of memory.");
  900.     exit(1);
  901.   }
  902.   /* Set the next block pointer to NULL. */
  903.   *(pool->firstblock) = (void *) NULL;
  904.   poolrestart(pool);
  905. }
  906.  
  907. /*****************************************************************************/
  908. /*                                                                           */
  909. /*  poolrestart()   Deallocate all items in a pool.                          */
  910. /*                                                                           */
  911. /*  The pool is returned to its starting state, except that no memory is     */
  912. /*  freed to the operating system.  Rather, the previously allocated blocks  */
  913. /*  are ready to be reused.                                                  */
  914. /*                                                                           */
  915. /*****************************************************************************/
  916.  
  917. void poolrestart(struct memorypool *pool)
  918. {
  919.   unsigned long alignptr;
  920.  
  921.   pool->items = 0;
  922.   pool->maxitems = 0;
  923.  
  924.   /* Set the currently active block. */
  925.   pool->nowblock = pool->firstblock;
  926.   /* Find the first item in the pool.  Increment by the size of (void *). */
  927.   alignptr = (unsigned long) (pool->nowblock + 1);
  928.   /* Align the item on an `alignbytes'-byte boundary. */
  929.   pool->nextitem = (void *)
  930.     (alignptr + (unsigned long) pool->alignbytes
  931.      - (alignptr % (unsigned long) pool->alignbytes));
  932.   /* There are lots of unallocated items left in this block. */
  933.   pool->unallocateditems = pool->itemsperblock;
  934.   /* The stack of deallocated items is empty. */
  935.   pool->deaditemstack = (void *) NULL;
  936. }
  937.  
  938. /*****************************************************************************/
  939. /*                                                                           */
  940. /*  pooldeinit()   Free to the operating system all memory taken by a pool.  */
  941. /*                                                                           */
  942. /*****************************************************************************/
  943.  
  944. void pooldeinit(struct memorypool *pool)
  945. {
  946.   while (pool->firstblock != (void **) NULL) {
  947.     pool->nowblock = (void **) *(pool->firstblock);
  948.     free(pool->firstblock);
  949.     pool->firstblock = pool->nowblock;
  950.   }
  951. }
  952.  
  953. /*****************************************************************************/
  954. /*                                                                           */
  955. /*  poolalloc()   Allocate space for an item.                                */
  956. /*                                                                           */
  957. /*****************************************************************************/
  958.  
  959. void *poolalloc(struct memorypool *pool)
  960. {
  961.   void *newitem;
  962.   void **newblock;
  963.   unsigned long alignptr;
  964.  
  965.   /* First check the linked list of dead items.  If the list is not   */
  966.   /*   empty, allocate an item from the list rather than a fresh one. */
  967.   if (pool->deaditemstack != (void *) NULL) {
  968.     newitem = pool->deaditemstack;               /* Take first item in list. */
  969.     pool->deaditemstack = * (void **) pool->deaditemstack;
  970.   } else {
  971.     /* Check if there are any free items left in the current block. */
  972.     if (pool->unallocateditems == 0) {
  973.       /* Check if another block must be allocated. */
  974.       if (*(pool->nowblock) == (void *) NULL) {
  975.         /* Allocate a new block of items, pointed to by the previous block. */
  976.         newblock = (void **) malloc(pool->itemsperblock * pool->itembytes
  977.                                     + sizeof(void *) + pool->alignbytes);
  978.         if (newblock == (void **) NULL) {
  979.           vTrace("Error:  Out of memory.");
  980.           exit(1);
  981.         }
  982.         *(pool->nowblock) = (void *) newblock;
  983.         /* The next block pointer is NULL. */
  984.         *newblock = (void *) NULL;
  985.       }
  986.       /* Move to the new block. */
  987.       pool->nowblock = (void **) *(pool->nowblock);
  988.       /* Find the first item in the block.    */
  989.       /*   Increment by the size of (void *). */
  990.       alignptr = (unsigned long) (pool->nowblock + 1);
  991.       /* Align the item on an `alignbytes'-byte boundary. */
  992.       pool->nextitem = (void *)
  993.         (alignptr + (unsigned long) pool->alignbytes
  994.          - (alignptr % (unsigned long) pool->alignbytes));
  995.       /* There are lots of unallocated items left in this block. */
  996.       pool->unallocateditems = pool->itemsperblock;
  997.     }
  998.     /* Allocate a new item. */
  999.     newitem = pool->nextitem;
  1000.     /* Advance `nextitem' pointer to next free item in block. */
  1001.     if (pool->itemwordtype == POINTER) {
  1002.       pool->nextitem = (void *) ((void **) pool->nextitem + pool->itemwords);
  1003.     } else {
  1004.       pool->nextitem = (void *) ((double *) pool->nextitem + pool->itemwords);
  1005.     }
  1006.     pool->unallocateditems--;
  1007.     pool->maxitems++;
  1008.   }
  1009.   pool->items++;
  1010.   return newitem;
  1011. }
  1012.  
  1013. /*****************************************************************************/
  1014. /*                                                                           */
  1015. /*  pooldealloc()   Deallocate space for an item.                            */
  1016. /*                                                                           */
  1017. /*  The deallocated space is stored in a queue for later reuse.              */
  1018. /*                                                                           */
  1019. /*****************************************************************************/
  1020.  
  1021. void pooldealloc(struct memorypool *pool, void* dyingitem)
  1022. {
  1023.   /* Push freshly killed item onto stack. */
  1024.   *((void **) dyingitem) = pool->deaditemstack;
  1025.   pool->deaditemstack = dyingitem;
  1026.   pool->items--;
  1027. }
  1028.  
  1029. /*****************************************************************************/
  1030. /*                                                                           */
  1031. /*  traversalinit()   Prepare to traverse the entire list of items.          */
  1032. /*                                                                           */
  1033. /*  This routine is used in conjunction with traverse().                     */
  1034. /*                                                                           */
  1035. /*****************************************************************************/
  1036.  
  1037. void traversalinit(struct memorypool *pool)
  1038. {
  1039.   unsigned long alignptr;
  1040.  
  1041.   /* Begin the traversal in the first block. */
  1042.   pool->pathblock = pool->firstblock;
  1043.   /* Find the first item in the block.  Increment by the size of (void *). */
  1044.   alignptr = (unsigned long) (pool->pathblock + 1);
  1045.   /* Align with item on an `alignbytes'-byte boundary. */
  1046.   pool->pathitem = (void *)
  1047.     (alignptr + (unsigned long) pool->alignbytes
  1048.      - (alignptr % (unsigned long) pool->alignbytes));
  1049.   /* Set the number of items left in the current block. */
  1050.   pool->pathitemsleft = pool->itemsperblock;
  1051. }
  1052.  
  1053. /*****************************************************************************/
  1054. /*                                                                           */
  1055. /*  traverse()   Find the next item in the list.                             */
  1056. /*                                                                           */
  1057. /*  This routine is used in conjunction with traversalinit().  Be forewarned */
  1058. /*  that this routine successively returns all items in the list, including  */
  1059. /*  deallocated ones on the deaditemqueue.  It's up to you to figure out     */
  1060. /*  which ones are actually dead.  Why?  I don't want to allocate extra      */
  1061. /*  space just to demarcate dead items.  It can usually be done more         */
  1062. /*  space-efficiently by a routine that knows something about the structure  */
  1063. /*  of the item.                                                             */
  1064. /*                                                                           */
  1065. /*****************************************************************************/
  1066.  
  1067. void *traverse(struct memorypool *pool)
  1068. {
  1069.   void *newitem;
  1070.   unsigned long alignptr;
  1071.  
  1072.   /* Stop upon exhausting the list of items. */
  1073.   if (pool->pathitem == pool->nextitem) {
  1074.     return (void *) NULL;
  1075.   }
  1076.   /* Check whether any untraversed items remain in the current block. */
  1077.   if (pool->pathitemsleft == 0) {
  1078.     /* Find the next block. */
  1079.     pool->pathblock = (void **) *(pool->pathblock);
  1080.     /* Find the first item in the block.  Increment by the size of (void *). */
  1081.     alignptr = (unsigned long) (pool->pathblock + 1);
  1082.     /* Align with item on an `alignbytes'-byte boundary. */
  1083.     pool->pathitem = (void *)
  1084.       (alignptr + (unsigned long) pool->alignbytes
  1085.        - (alignptr % (unsigned long) pool->alignbytes));
  1086.     /* Set the number of items left in the current block. */
  1087.     pool->pathitemsleft = pool->itemsperblock;
  1088.   }
  1089.   newitem = pool->pathitem;
  1090.   /* Find the next item in the block. */
  1091.   if (pool->itemwordtype == POINTER) {
  1092.     pool->pathitem = (void *) ((void **) pool->pathitem + pool->itemwords);
  1093.   } else {
  1094.     pool->pathitem = (void *) ((double *) pool->pathitem + pool->itemwords);
  1095.   }
  1096.   pool->pathitemsleft--;
  1097.   return newitem;
  1098. }
  1099.  
  1100. /*****************************************************************************/
  1101. /*                                                                           */
  1102. /*  dummyinit()   Initialize the triangle that fills "outer space" and the   */
  1103. /*                omnipresent shell edge.                                    */
  1104. /*                                                                           */
  1105. /*  The triangle that fills "outer space", called `dummytri', is pointed to  */
  1106. /*  by every triangle and shell edge on a boundary (be it outer or inner) of */
  1107. /*  the triangulation.  Also, `dummytri' points to one of the triangles on   */
  1108. /*  the convex hull (until the holes and concavities are carved), making it  */
  1109. /*  possible to find a starting triangle for point location.                 */
  1110. /*                                                                           */
  1111. /*  The omnipresent shell edge, `dummysh', is pointed to by every triangle   */
  1112. /*  or shell edge that doesn't have a full complement of real shell edges    */
  1113. /*  to point to.                                                             */
  1114. /*                                                                           */
  1115. /*****************************************************************************/
  1116.  
  1117. void dummyinit(int trianglewords, int shellewords)
  1118. {
  1119.   unsigned long alignptr;
  1120.  
  1121.   /* `triwords' and `shwords' are used by the mesh manipulation primitives */
  1122.   /*   to extract orientations of triangles and shell edges from pointers. */
  1123.   triwords = trianglewords;       /* Initialize `triwords' once and for all. */
  1124.   shwords = shellewords;           /* Initialize `shwords' once and for all. */
  1125.  
  1126.   /* Set up `dummytri', the `triangle' that occupies "outer space". */
  1127.   dummytribase = (triangle *) malloc(triwords * sizeof(triangle)
  1128.                                      + triangles.alignbytes);
  1129.   if (dummytribase == (triangle *) NULL) {
  1130.     vTrace("Error:  Out of memory.");
  1131.     exit(1);
  1132.   }
  1133.   /* Align `dummytri' on a `triangles.alignbytes'-byte boundary. */
  1134.   alignptr = (unsigned long) dummytribase;
  1135.   dummytri = (triangle *)
  1136.     (alignptr + (unsigned long) triangles.alignbytes
  1137.      - (alignptr % (unsigned long) triangles.alignbytes));
  1138.   /* Initialize the three adjoining triangles to be "outer space".  These  */
  1139.   /*   will eventually be changed by various bonding operations, but their */
  1140.   /*   values don't really matter, as long as they can legally be          */
  1141.   /*   dereferenced.                                                       */
  1142.   dummytri[0] = (triangle) dummytri;
  1143.   dummytri[1] = (triangle) dummytri;
  1144.   dummytri[2] = (triangle) dummytri;
  1145.   /* Three NULL vertex points. */
  1146.   dummytri[3] = (triangle) NULL;
  1147.   dummytri[4] = (triangle) NULL;
  1148.   dummytri[5] = (triangle) NULL;
  1149.  
  1150.   if (useshelles) {
  1151.     /* Set up `dummysh', the omnipresent "shell edge" pointed to by any      */
  1152.     /*   triangle side or shell edge end that isn't attached to a real shell */
  1153.     /*   edge.                                                               */
  1154.     dummyshbase = (shelle *) malloc(shwords * sizeof(shelle)
  1155.                                     + shelles.alignbytes);
  1156.     if (dummyshbase == (shelle *) NULL) {
  1157.       vTrace("Error:  Out of memory.");
  1158.       exit(1);
  1159.     }
  1160.     /* Align `dummysh' on a `shelles.alignbytes'-byte boundary. */
  1161.     alignptr = (unsigned long) dummyshbase;
  1162.     dummysh = (shelle *)
  1163.       (alignptr + (unsigned long) shelles.alignbytes
  1164.        - (alignptr % (unsigned long) shelles.alignbytes));
  1165.     /* Initialize the two adjoining shell edges to be the omnipresent shell */
  1166.     /*   edge.  These will eventually be changed by various bonding         */
  1167.     /*   operations, but their values don't really matter, as long as they  */
  1168.     /*   can legally be dereferenced.                                       */
  1169.     dummysh[0] = (shelle) dummysh;
  1170.     dummysh[1] = (shelle) dummysh;
  1171.     /* Two NULL vertex points. */
  1172.     dummysh[2] = (shelle) NULL;
  1173.     dummysh[3] = (shelle) NULL;
  1174.     /* Initialize the two adjoining triangles to be "outer space". */
  1175.     dummysh[4] = (shelle) dummytri;
  1176.     dummysh[5] = (shelle) dummytri;
  1177.     /* Set the boundary marker to zero. */
  1178.     * (int *) (dummysh + 6) = 0;
  1179.  
  1180.     /* Initialize the three adjoining shell edges of `dummytri' to be */
  1181.     /*   the omnipresent shell edge.                                  */
  1182.     dummytri[6] = (triangle) dummysh;
  1183.     dummytri[7] = (triangle) dummysh;
  1184.     dummytri[8] = (triangle) dummysh;
  1185.   }
  1186. }
  1187.  
  1188. /*****************************************************************************/
  1189. /*                                                                           */
  1190. /*  initializepointpool()   Calculate the size of the point data structure   */
  1191. /*                          and initialize its memory pool.                  */
  1192. /*                                                                           */
  1193. /*  This routine also computes the `pointmarkindex' and `point2triindex'     */
  1194. /*  indices used to find values within each point.                           */
  1195. /*                                                                           */
  1196. /*****************************************************************************/
  1197.  
  1198. void initializepointpool(void)
  1199. {
  1200.   int pointsize;
  1201.  
  1202.   /* The index within each point at which the boundary marker is found.  */
  1203.   /*   Ensure the point marker is aligned to a sizeof(int)-byte address. */
  1204.   pointmarkindex = ((mesh_dim + nextras) * sizeof(double) + sizeof(int) - 1)
  1205.                  / sizeof(int);
  1206.   pointsize = (pointmarkindex + 1) * sizeof(int);
  1207.     /* The index within each point at which a triangle pointer is found.   */
  1208.     /*   Ensure the pointer is aligned to a sizeof(triangle)-byte address. */
  1209.     point2triindex = (pointsize + sizeof(triangle) - 1) / sizeof(triangle);
  1210.     pointsize = (point2triindex + 1) * sizeof(triangle);
  1211.  
  1212.     /* Initialize the pool of points. */
  1213.   poolinit(&points, pointsize, POINTPERBLOCK,
  1214.            (sizeof(double) >= sizeof(triangle)) ? FLOATINGPOINT : POINTER, 0);
  1215. }
  1216.  
  1217. /*****************************************************************************/
  1218. /*                                                                           */
  1219. /*  initializetrisegpools()   Calculate the sizes of the triangle and shell  */
  1220. /*                            edge data structures and initialize their      */
  1221. /*                            memory pools.                                  */
  1222. /*                                                                           */
  1223. /*  This routine also computes the `highorderindex', `elemattribindex', and  */
  1224. /*  `areaboundindex' indices used to find values within each triangle.       */
  1225. /*                                                                           */
  1226. /*****************************************************************************/
  1227.  
  1228. void initializetrisegpools(void)
  1229. {
  1230.   int trisize;
  1231.  
  1232.   /* The index within each triangle at which the extra nodes (above three)  */
  1233.   /*   associated with high order elements are found.  There are three      */
  1234.   /*   pointers to other triangles, three pointers to corners, and possibly */
  1235.   /*   three pointers to shell edges before the extra nodes.                */
  1236.   highorderindex = 6 + (useshelles * 3);
  1237.   /* The number of bytes occupied by a triangle. */
  1238.   trisize = ((order + 1) * (order + 2) / 2 + (highorderindex - 3)) *
  1239.             sizeof(triangle);
  1240.   /* The index within each triangle at which its attributes are found, */
  1241.   /*   where the index is measured in REALs.                           */
  1242.   elemattribindex = (trisize + sizeof(double) - 1) / sizeof(double);
  1243.   /* The index within each triangle at which the maximum area constraint  */
  1244.   /*   is found, where the index is measured in REALs.*/
  1245.   areaboundindex = elemattribindex + eextras;
  1246.   /* If triangle attributes or an area bound are needed, increase the number */
  1247.   /*   of bytes occupied by a triangle.                                      */
  1248.  if (eextras > 0) {
  1249.     trisize = areaboundindex * sizeof(double);
  1250.   }
  1251.  
  1252.   /* Having determined the memory size of a triangle, initialize the pool. */
  1253.   poolinit(&triangles, trisize, TRIPERBLOCK, POINTER, 4);
  1254.  
  1255.   if (useshelles) {
  1256.     /* Initialize the pool of shell edges. */
  1257.     poolinit(&shelles, 6 * sizeof(triangle) + sizeof(int), SHELLEPERBLOCK,
  1258.              POINTER, 4);
  1259.  
  1260.     /* Initialize the "outer space" triangle and omnipresent shell edge. */
  1261.     dummyinit(triangles.itemwords, shelles.itemwords);
  1262.   } else {
  1263.     /* Initialize the "outer space" triangle. */
  1264.     dummyinit(triangles.itemwords, 0);
  1265.   }
  1266. }
  1267.  
  1268. /*****************************************************************************/
  1269. /*                                                                           */
  1270. /*  triangledealloc()   Deallocate space for a triangle, marking it dead.    */
  1271. /*                                                                           */
  1272. /*****************************************************************************/
  1273.  
  1274. void triangledealloc(triangle *dyingtriangle)
  1275. {
  1276.   /* Set triangle's vertices to NULL.  This makes it possible to        */
  1277.   /*   detect dead triangles when traversing the list of all triangles. */
  1278.   dyingtriangle[3] = (triangle) NULL;
  1279.   dyingtriangle[4] = (triangle) NULL;
  1280.   dyingtriangle[5] = (triangle) NULL;
  1281.   pooldealloc(&triangles, (void *) dyingtriangle);
  1282. }
  1283.  
  1284. /*****************************************************************************/
  1285. /*                                                                           */
  1286. /*  triangletraverse()   Traverse the triangles, skipping dead ones.         */
  1287. /*                                                                           */
  1288. /*****************************************************************************/
  1289.  
  1290. triangle *triangletraverse(void)
  1291. {
  1292.   triangle *newtriangle;
  1293.  
  1294.   do {
  1295.     newtriangle = (triangle *) traverse(&triangles);
  1296.     if (newtriangle == (triangle *) NULL) {
  1297.       return (triangle *) NULL;
  1298.     }
  1299.   } while (newtriangle[3] == (triangle) NULL);            /* Skip dead ones. */
  1300.   return newtriangle;
  1301. }
  1302.  
  1303. /*****************************************************************************/
  1304. /*                                                                           */
  1305. /*  shelledealloc()   Deallocate space for a shell edge, marking it dead.    */
  1306. /*                                                                           */
  1307. /*****************************************************************************/
  1308.  
  1309. void shelledealloc(shelle *dyingshelle)
  1310. {
  1311.   /* Set shell edge's vertices to NULL.  This makes it possible to */
  1312.   /*   detect dead shells when traversing the list of all shells.  */
  1313.   dyingshelle[2] = (shelle) NULL;
  1314.   dyingshelle[3] = (shelle) NULL;
  1315.   pooldealloc(&shelles, (void *) dyingshelle);
  1316. }
  1317.  
  1318. /*****************************************************************************/
  1319. /*                                                                           */
  1320. /*  pointdealloc()   Deallocate space for a point, marking it dead.          */
  1321. /*                                                                           */
  1322. /*****************************************************************************/
  1323.  
  1324. void pointdealloc(point dyingpoint)
  1325. {
  1326.   /* Mark the point as dead.  This makes it possible to detect dead points */
  1327.   /*   when traversing the list of all points.                             */
  1328.   setpointmark(dyingpoint, DEADPOINT);
  1329.   pooldealloc(&points, (void *) dyingpoint);
  1330. }
  1331.  
  1332. /*****************************************************************************/
  1333. /*                                                                           */
  1334. /*  pointtraverse()   Traverse the points, skipping dead ones.               */
  1335. /*                                                                           */
  1336. /*****************************************************************************/
  1337.  
  1338. point pointtraverse(void)
  1339. {
  1340.   point newpoint;
  1341.  
  1342.   do {
  1343.     newpoint = (point) traverse(&points);
  1344.     if (newpoint == (point) NULL) {
  1345.       return (point) NULL;
  1346.     }
  1347.   } while (pointmark(newpoint) == DEADPOINT);             /* Skip dead ones. */
  1348.   return newpoint;
  1349. }
  1350.  
  1351. /*****************************************************************************/
  1352. /*                                                                           */
  1353. /*  getpoint()   Get a specific point, by number, from the list.             */
  1354. /*                                                                           */
  1355. /*  The first point is number 0.                                 */
  1356. /*                                                                           */
  1357. /*  Note that this takes O(n) time (with a small constant, if POINTPERBLOCK  */
  1358. /*  is large).  I don't care to take the trouble to make it work in constant */
  1359. /*  time.                                                                    */
  1360. /*                                                                           */
  1361. /*****************************************************************************/
  1362.  
  1363. point getpoint(int number)
  1364. {
  1365.   void **getblock;
  1366.   point foundpoint;
  1367.   unsigned long alignptr;
  1368.   int current;
  1369.  
  1370.   getblock = points.firstblock;
  1371.   current = 0;
  1372.   /* Find the right block. */
  1373.   while (current + points.itemsperblock <= number) {
  1374.     getblock = (void **) *getblock;
  1375.     current += points.itemsperblock;
  1376.   }
  1377.   /* Now find the right point. */
  1378.   alignptr = (unsigned long) (getblock + 1);
  1379.   foundpoint = (point) (alignptr + (unsigned long) points.alignbytes
  1380.                         - (alignptr % (unsigned long) points.alignbytes));
  1381.   while (current < number) {
  1382.     foundpoint += points.itemwords;
  1383.     current++;
  1384.   }
  1385.   return foundpoint;
  1386. }
  1387.  
  1388. /*****************************************************************************/
  1389. /*                                                                           */
  1390. /*  triangledeinit()   Free all remaining allocated memory.                  */
  1391. /*                                                                           */
  1392. /*****************************************************************************/
  1393.  
  1394. void triangledeinit(void)
  1395. {
  1396.   pooldeinit(&triangles);
  1397.   free(dummytribase);
  1398.   if (useshelles) {
  1399.     pooldeinit(&shelles);
  1400.     free(dummyshbase);
  1401.   }
  1402.   pooldeinit(&points);
  1403. }
  1404.  
  1405. /**                                                                         **/
  1406. /**                                                                         **/
  1407. /********* Memory management routines end here                       *********/
  1408.  
  1409. /********* Constructors begin here                                   *********/
  1410. /**                                                                         **/
  1411. /**                                                                         **/
  1412.  
  1413. /*****************************************************************************/
  1414. /*                                                                           */
  1415. /*  maketriangle()   Create a new triangle with orientation zero.            */
  1416. /*                                                                           */
  1417. /*****************************************************************************/
  1418.  
  1419. void maketriangle(struct triedge *newtriedge)
  1420. {
  1421.   int i;
  1422.  
  1423.   newtriedge->tri = (triangle *) poolalloc(&triangles);
  1424.   /* Initialize the three adjoining triangles to be "outer space". */
  1425.   newtriedge->tri[0] = (triangle) dummytri;
  1426.   newtriedge->tri[1] = (triangle) dummytri;
  1427.   newtriedge->tri[2] = (triangle) dummytri;
  1428.   /* Three NULL vertex points. */
  1429.   newtriedge->tri[3] = (triangle) NULL;
  1430.   newtriedge->tri[4] = (triangle) NULL;
  1431.   newtriedge->tri[5] = (triangle) NULL;
  1432.   /* Initialize the three adjoining shell edges to be the omnipresent */
  1433.   /*   shell edge.                                                    */
  1434.   if (useshelles) {
  1435.     newtriedge->tri[6] = (triangle) dummysh;
  1436.     newtriedge->tri[7] = (triangle) dummysh;
  1437.     newtriedge->tri[8] = (triangle) dummysh;
  1438.   }
  1439.   for (i = 0; i < eextras; i++) {
  1440.     setelemattribute(*newtriedge, i, 0.0);
  1441.   }
  1442.  
  1443.   newtriedge->orient = 0;
  1444. }
  1445.  
  1446. /*****************************************************************************/
  1447. /*                                                                           */
  1448. /*  makeshelle()   Create a new shell edge with orientation zero.            */
  1449. /*                                                                           */
  1450. /*****************************************************************************/
  1451.  
  1452. void makeshelle(struct edge *newedge)
  1453. {
  1454.   newedge->sh = (shelle *) poolalloc(&shelles);
  1455.   /* Initialize the two adjoining shell edges to be the omnipresent */
  1456.   /*   shell edge.                                                  */
  1457.   newedge->sh[0] = (shelle) dummysh;
  1458.   newedge->sh[1] = (shelle) dummysh;
  1459.   /* Two NULL vertex points. */
  1460.   newedge->sh[2] = (shelle) NULL;
  1461.   newedge->sh[3] = (shelle) NULL;
  1462.   /* Initialize the two adjoining triangles to be "outer space". */
  1463.   newedge->sh[4] = (shelle) dummytri;
  1464.   newedge->sh[5] = (shelle) dummytri;
  1465.   /* Set the boundary marker to zero. */
  1466.   setmark(*newedge, 0);
  1467.  
  1468.   newedge->shorient = 0;
  1469. }
  1470.  
  1471. /**                                                                         **/
  1472. /**                                                                         **/
  1473. /********* Constructors end here                                     *********/
  1474.  
  1475. /********* Determinant evaluation routines begin here                *********/
  1476. /**                                                                         **/
  1477. /**                                                                         **/
  1478.  
  1479. /* The adaptive exact arithmetic geometric predicates implemented herein are */
  1480. /*   described in detail in my Technical Report CMU-CS-96-140.  The complete */
  1481. /*   reference is given in the header.                                       */
  1482.  
  1483. /* Which of the following two methods of finding the absolute values is      */
  1484. /*   fastest is compiler-dependent.  A few compilers can inline and optimize */
  1485. /*   the fabs() call; but most will incur the overhead of a function call,   */
  1486. /*   which is disastrously slow.  A faster way on IEEE machines might be to  */
  1487. /*   mask the appropriate bit, but that's difficult to do in C.              */
  1488.  
  1489. // #define Absolute(a)  ((a) >= 0.0 ? (a) : -(a))
  1490. #define Absolute(a)  fabs(a)
  1491.  
  1492. /* Many of the operations are broken up into two pieces, a main part that    */
  1493. /*   performs an approximate operation, and a "tail" that computes the       */
  1494. /*   roundoff error of that operation.                                       */
  1495. /*                                                                           */
  1496. /* The operations Fast_Two_Sum(), Fast_Two_Diff(), Two_Sum(), Two_Diff(),    */
  1497. /*   Split(), and Two_Product() are all implemented as described in the      */
  1498. /*   reference.  Each of these macros requires certain variables to be       */
  1499. /*   defined in the calling routine.  The variables `bvirt', `c', `abig',    */
  1500. /*   `_i', `_j', `_k', `_l', `_m', and `_n' are declared `INEXACT' because   */
  1501. /*   they store the result of an operation that may incur roundoff error.    */
  1502. /*   The input parameter `x' (or the highest numbered `x_' parameter) must   */
  1503. /*   also be declared `INEXACT'.                                             */
  1504.  
  1505. #define Fast_Two_Sum_Tail(a, b, x, y) \
  1506.   bvirt = x - a; \
  1507.   y = b - bvirt
  1508.  
  1509. #define Fast_Two_Sum(a, b, x, y) \
  1510.   x = (double) (a + b); \
  1511.   Fast_Two_Sum_Tail(a, b, x, y)
  1512.  
  1513. #define Two_Sum_Tail(a, b, x, y) \
  1514.   bvirt = (double) (x - a); \
  1515.   avirt = x - bvirt; \
  1516.   bround = b - bvirt; \
  1517.   around = a - avirt; \
  1518.   y = around + bround
  1519.  
  1520. #define Two_Sum(a, b, x, y) \
  1521.   x = (double) (a + b); \
  1522.   Two_Sum_Tail(a, b, x, y)
  1523.  
  1524. #define Two_Diff_Tail(a, b, x, y) \
  1525.   bvirt = (double) (a - x); \
  1526.   avirt = x + bvirt; \
  1527.   bround = bvirt - b; \
  1528.   around = a - avirt; \
  1529.   y = around + bround
  1530.  
  1531. #define Two_Diff(a, b, x, y) \
  1532.   x = (double) (a - b); \
  1533.   Two_Diff_Tail(a, b, x, y)
  1534.  
  1535. #define Split(a, ahi, alo) \
  1536.   c = (double) (splitter * a); \
  1537.   abig = (double) (c - a); \
  1538.   ahi = c - abig; \
  1539.   alo = a - ahi
  1540.  
  1541. #define Two_Product_Tail(a, b, x, y) \
  1542.   Split(a, ahi, alo); \
  1543.   Split(b, bhi, blo); \
  1544.   err1 = x - (ahi * bhi); \
  1545.   err2 = err1 - (alo * bhi); \
  1546.   err3 = err2 - (ahi * blo); \
  1547.   y = (alo * blo) - err3
  1548.  
  1549. #define Two_Product(a, b, x, y) \
  1550.   x = (double) (a * b); \
  1551.   Two_Product_Tail(a, b, x, y)
  1552.  
  1553. /* Two_Product_Presplit() is Two_Product() where one of the inputs has       */
  1554. /*   already been split.  Avoids redundant splitting.                        */
  1555.  
  1556. #define Two_Product_Presplit(a, b, bhi, blo, x, y) \
  1557.   x = (double) (a * b); \
  1558.   Split(a, ahi, alo); \
  1559.   err1 = x - (ahi * bhi); \
  1560.   err2 = err1 - (alo * bhi); \
  1561.   err3 = err2 - (ahi * blo); \
  1562.   y = (alo * blo) - err3
  1563.  
  1564. /* Square() can be done more quickly than Two_Product().                     */
  1565.  
  1566. #define Square_Tail(a, x, y) \
  1567.   Split(a, ahi, alo); \
  1568.   err1 = x - (ahi * ahi); \
  1569.   err3 = err1 - ((ahi + ahi) * alo); \
  1570.   y = (alo * alo) - err3
  1571.  
  1572. #define Square(a, x, y) \
  1573.   x = (double) (a * a); \
  1574.   Square_Tail(a, x, y)
  1575.  
  1576. /* Macros for summing expansions of various fixed lengths.  These are all    */
  1577. /*   unrolled versions of Expansion_Sum().                                   */
  1578.  
  1579. #define Two_One_Sum(a1, a0, b, x2, x1, x0) \
  1580.   Two_Sum(a0, b , _i, x0); \
  1581.   Two_Sum(a1, _i, x2, x1)
  1582.  
  1583. #define Two_One_Diff(a1, a0, b, x2, x1, x0) \
  1584.   Two_Diff(a0, b , _i, x0); \
  1585.   Two_Sum( a1, _i, x2, x1)
  1586.  
  1587. #define Two_Two_Sum(a1, a0, b1, b0, x3, x2, x1, x0) \
  1588.   Two_One_Sum(a1, a0, b0, _j, _0, x0); \
  1589.   Two_One_Sum(_j, _0, b1, x3, x2, x1)
  1590.  
  1591. #define Two_Two_Diff(a1, a0, b1, b0, x3, x2, x1, x0) \
  1592.   Two_One_Diff(a1, a0, b0, _j, _0, x0); \
  1593.   Two_One_Diff(_j, _0, b1, x3, x2, x1)
  1594.  
  1595. /*****************************************************************************/
  1596. /*                                                                           */
  1597. /*  exactinit()   Initialize the variables used for exact arithmetic.        */
  1598. /*                                                                           */
  1599. /*  `epsilon' is the largest power of two such that 1.0 + epsilon = 1.0 in   */
  1600. /*  floating-point arithmetic.  `epsilon' bounds the relative roundoff       */
  1601. /*  error.  It is used for floating-point error analysis.                    */
  1602. /*                                                                           */
  1603. /*  `splitter' is used to split floating-point numbers into two half-        */
  1604. /*  length significands for exact multiplication.                            */
  1605. /*                                                                           */
  1606. /*  I imagine that a highly optimizing compiler might be too smart for its   */
  1607. /*  own good, and somehow cause this routine to fail, if it pretends that    */
  1608. /*  floating-point arithmetic is too much like real arithmetic.              */
  1609. /*                                                                           */
  1610. /*  Don't change this routine unless you fully understand it.                */
  1611. /*                                                                           */
  1612. /*****************************************************************************/
  1613.  
  1614. void exactinit(void)
  1615. {
  1616.   double half;
  1617.   double check, lastcheck;
  1618.   int every_other;
  1619.  
  1620.   every_other = 1;
  1621.   half = 0.5;
  1622.   epsilon = 1.0;
  1623.   splitter = 1.0;
  1624.   check = 1.0;
  1625.   /* Repeatedly divide `epsilon' by two until it is too small to add to      */
  1626.   /*   one without causing roundoff.  (Also check if the sum is equal to     */
  1627.   /*   the previous sum, for machines that round up instead of using exact   */
  1628.   /*   rounding.  Not that these routines will work on such machines anyway. */
  1629.   do {
  1630.     lastcheck = check;
  1631.     epsilon *= half;
  1632.     if (every_other) {
  1633.       splitter *= 2.0;
  1634.     }
  1635.     every_other = !every_other;
  1636.     check = 1.0 + epsilon;
  1637.   } while ((check != 1.0) && (check != lastcheck));
  1638.   splitter += 1.0;
  1639.  
  1640.   /* Error bounds for orientation and incircle tests. */
  1641.   resulterrbound = (3.0 + 8.0 * epsilon) * epsilon;
  1642.   ccwerrboundA = (3.0 + 16.0 * epsilon) * epsilon;
  1643.   ccwerrboundB = (2.0 + 12.0 * epsilon) * epsilon;
  1644.   ccwerrboundC = (9.0 + 64.0 * epsilon) * epsilon * epsilon;
  1645.   iccerrboundA = (10.0 + 96.0 * epsilon) * epsilon;
  1646.   iccerrboundB = (4.0 + 48.0 * epsilon) * epsilon;
  1647.   iccerrboundC = (44.0 + 576.0 * epsilon) * epsilon * epsilon;
  1648. }
  1649.  
  1650. /*****************************************************************************/
  1651. /*                                                                           */
  1652. /*  fast_expansion_sum_zeroelim()   Sum two expansions, eliminating zero     */
  1653. /*                                  components from the output expansion.    */
  1654. /*                                                                           */
  1655. /*  Sets h = e + f.  See my Robust Predicates paper for details.             */
  1656. /*                                                                           */
  1657. /*  If round-to-even is used (as with IEEE 754), maintains the strongly      */
  1658. /*  nonoverlapping property.  (That is, if e is strongly nonoverlapping, h   */
  1659. /*  will be also.)  Does NOT maintain the nonoverlapping or nonadjacent      */
  1660. /*  properties.                                                              */
  1661. /*                                                                           */
  1662. /*****************************************************************************/
  1663.  
  1664. int fast_expansion_sum_zeroelim(int elen,double *e,int flen,double *f,double *h)  /* h cannot be e or f. */
  1665. {
  1666.   double Q;
  1667.   INEXACT double Qnew;
  1668.   INEXACT double hh;
  1669.   INEXACT double bvirt;
  1670.   double avirt, bround, around;
  1671.   int eindex, findex, hindex;
  1672.   double enow, fnow;
  1673.  
  1674.   enow = e[0];
  1675.   fnow = f[0];
  1676.   eindex = findex = 0;
  1677.   if ((fnow > enow) == (fnow > -enow)) {
  1678.     Q = enow;
  1679.     enow = e[++eindex];
  1680.   } else {
  1681.     Q = fnow;
  1682.     fnow = f[++findex];
  1683.   }
  1684.   hindex = 0;
  1685.   if ((eindex < elen) && (findex < flen)) {
  1686.     if ((fnow > enow) == (fnow > -enow)) {
  1687.       Fast_Two_Sum(enow, Q, Qnew, hh);
  1688.       enow = e[++eindex];
  1689.     } else {
  1690.       Fast_Two_Sum(fnow, Q, Qnew, hh);
  1691.       fnow = f[++findex];
  1692.     }
  1693.     Q = Qnew;
  1694.     if (hh != 0.0) {
  1695.       h[hindex++] = hh;
  1696.     }
  1697.     while ((eindex < elen) && (findex < flen)) {
  1698.       if ((fnow > enow) == (fnow > -enow)) {
  1699.         Two_Sum(Q, enow, Qnew, hh);
  1700.         enow = e[++eindex];
  1701.       } else {
  1702.         Two_Sum(Q, fnow, Qnew, hh);
  1703.         fnow = f[++findex];
  1704.       }
  1705.       Q = Qnew;
  1706.       if (hh != 0.0) {
  1707.         h[hindex++] = hh;
  1708.       }
  1709.     }
  1710.   }
  1711.   while (eindex < elen) {
  1712.     Two_Sum(Q, enow, Qnew, hh);
  1713.     enow = e[++eindex];
  1714.     Q = Qnew;
  1715.     if (hh != 0.0) {
  1716.       h[hindex++] = hh;
  1717.     }
  1718.   }
  1719.   while (findex < flen) {
  1720.     Two_Sum(Q, fnow, Qnew, hh);
  1721.     fnow = f[++findex];
  1722.     Q = Qnew;
  1723.     if (hh != 0.0) {
  1724.       h[hindex++] = hh;
  1725.     }
  1726.   }
  1727.   if ((Q != 0.0) || (hindex == 0)) {
  1728.     h[hindex++] = Q;
  1729.   }
  1730.   return hindex;
  1731. }
  1732.  
  1733. /*****************************************************************************/
  1734. /*                                                                           */
  1735. /*  scale_expansion_zeroelim()   Multiply an expansion by a scalar,          */
  1736. /*                               eliminating zero components from the        */
  1737. /*                               output expansion.                           */
  1738. /*                                                                           */
  1739. /*  Sets h = be.  See my Robust Predicates paper for details.                */
  1740. /*                                                                           */
  1741. /*  Maintains the nonoverlapping property.  If round-to-even is used (as     */
  1742. /*  with IEEE 754), maintains the strongly nonoverlapping and nonadjacent    */
  1743. /*  properties as well.  (That is, if e has one of these properties, so      */
  1744. /*  will h.)                                                                 */
  1745. /*                                                                           */
  1746. /*****************************************************************************/
  1747.  
  1748. int scale_expansion_zeroelim(int elen,double *e,double b,double *h)   /* e and h cannot be the same. */
  1749. {
  1750.   INEXACT double Q, sum;
  1751.   double hh;
  1752.   INEXACT double product1;
  1753.   double product0;
  1754.   int eindex, hindex;
  1755.   double enow;
  1756.   INEXACT double bvirt;
  1757.   double avirt, bround, around;
  1758.   INEXACT double c;
  1759.   INEXACT double abig;
  1760.   double ahi, alo, bhi, blo;
  1761.   double err1, err2, err3;
  1762.  
  1763.   Split(b, bhi, blo);
  1764.   Two_Product_Presplit(e[0], b, bhi, blo, Q, hh);
  1765.   hindex = 0;
  1766.   if (hh != 0) {
  1767.     h[hindex++] = hh;
  1768.   }
  1769.   for (eindex = 1; eindex < elen; eindex++) {
  1770.     enow = e[eindex];
  1771.     Two_Product_Presplit(enow, b, bhi, blo, product1, product0);
  1772.     Two_Sum(Q, product0, sum, hh);
  1773.     if (hh != 0) {
  1774.       h[hindex++] = hh;
  1775.     }
  1776.     Fast_Two_Sum(product1, sum, Q, hh);
  1777.     if (hh != 0) {
  1778.       h[hindex++] = hh;
  1779.     }
  1780.   }
  1781.   if ((Q != 0.0) || (hindex == 0)) {
  1782.     h[hindex++] = Q;
  1783.   }
  1784.   return hindex;
  1785. }
  1786.  
  1787. /*****************************************************************************/
  1788. /*                                                                           */
  1789. /*  estimate()   Produce a one-word estimate of an expansion's value.        */
  1790. /*                                                                           */
  1791. /*  See my Robust Predicates paper for details.                              */
  1792. /*                                                                           */
  1793. /*****************************************************************************/
  1794.  
  1795. double estimate(int elen, double *e)
  1796. {
  1797.   double Q;
  1798.   int eindex;
  1799.  
  1800.   Q = e[0];
  1801.   for (eindex = 1; eindex < elen; eindex++) {
  1802.     Q += e[eindex];
  1803.   }
  1804.   return Q;
  1805. }
  1806.  
  1807. /*****************************************************************************/
  1808. /*                                                                           */
  1809. /*  counterclockwise()   Return a positive value if the points pa, pb, and   */
  1810. /*                       pc occur in counterclockwise order; a negative      */
  1811. /*                       value if they occur in clockwise order; and zero    */
  1812. /*                       if they are collinear.  The result is also a rough  */
  1813. /*                       approximation of twice the signed area of the       */
  1814. /*                       triangle defined by the three points.               */
  1815. /*                                                                           */
  1816. /*  Uses exact arithmetic if necessary to ensure a correct answer.  The      */
  1817. /*  result returned is the determinant of a matrix.  This determinant is     */
  1818. /*  computed adaptively, in the sense that exact arithmetic is used only to  */
  1819. /*  the degree it is needed to ensure that the returned value has the        */
  1820. /*  correct sign.  Hence, this function is usually quite fast, but will run  */
  1821. /*  more slowly when the input points are collinear or nearly so.            */
  1822. /*                                                                           */
  1823. /*  See my Robust Predicates paper for details.                              */
  1824. /*                                                                           */
  1825. /*****************************************************************************/
  1826.  
  1827. double counterclockwiseadapt(
  1828. point pa,
  1829. point pb,
  1830. point pc,
  1831. double detsum)
  1832. {
  1833.   INEXACT double acx, acy, bcx, bcy;
  1834.   double acxtail, acytail, bcxtail, bcytail;
  1835.   INEXACT double detleft, detright;
  1836.   double detlefttail, detrighttail;
  1837.   double det, errbound;
  1838.   double B[4], C1[8], C2[12], D[16];
  1839.   INEXACT double B3;
  1840.   int C1length, C2length, Dlength;
  1841.   double u[4];
  1842.   INEXACT double u3;
  1843.   INEXACT double s1, t1;
  1844.   double s0, t0;
  1845.  
  1846.   INEXACT double bvirt;
  1847.   double avirt, bround, around;
  1848.   INEXACT double c;
  1849.   INEXACT double abig;
  1850.   double ahi, alo, bhi, blo;
  1851.   double err1, err2, err3;
  1852.   INEXACT double _i, _j;
  1853.   double _0;
  1854.  
  1855.   acx = (double) (pa[0] - pc[0]);
  1856.   bcx = (double) (pb[0] - pc[0]);
  1857.   acy = (double) (pa[1] - pc[1]);
  1858.   bcy = (double) (pb[1] - pc[1]);
  1859.  
  1860.   Two_Product(acx, bcy, detleft, detlefttail);
  1861.   Two_Product(acy, bcx, detright, detrighttail);
  1862.  
  1863.   Two_Two_Diff(detleft, detlefttail, detright, detrighttail,
  1864.                B3, B[2], B[1], B[0]);
  1865.   B[3] = B3;
  1866.  
  1867.   det = estimate(4, B);
  1868.   errbound = ccwerrboundB * detsum;
  1869.   if ((det >= errbound) || (-det >= errbound)) {
  1870.     return det;
  1871.   }
  1872.  
  1873.   Two_Diff_Tail(pa[0], pc[0], acx, acxtail);
  1874.   Two_Diff_Tail(pb[0], pc[0], bcx, bcxtail);
  1875.   Two_Diff_Tail(pa[1], pc[1], acy, acytail);
  1876.   Two_Diff_Tail(pb[1], pc[1], bcy, bcytail);
  1877.  
  1878.   if ((acxtail == 0.0) && (acytail == 0.0)
  1879.       && (bcxtail == 0.0) && (bcytail == 0.0)) {
  1880.     return det;
  1881.   }
  1882.  
  1883.   errbound = ccwerrboundC * detsum + resulterrbound * Absolute(det);
  1884.   det += (acx * bcytail + bcy * acxtail)
  1885.        - (acy * bcxtail + bcx * acytail);
  1886.   if ((det >= errbound) || (-det >= errbound)) {
  1887.     return det;
  1888.   }
  1889.  
  1890.   Two_Product(acxtail, bcy, s1, s0);
  1891.   Two_Product(acytail, bcx, t1, t0);
  1892.   Two_Two_Diff(s1, s0, t1, t0, u3, u[2], u[1], u[0]);
  1893.   u[3] = u3;
  1894.   C1length = fast_expansion_sum_zeroelim(4, B, 4, u, C1);
  1895.  
  1896.   Two_Product(acx, bcytail, s1, s0);
  1897.   Two_Product(acy, bcxtail, t1, t0);
  1898.   Two_Two_Diff(s1, s0, t1, t0, u3, u[2], u[1], u[0]);
  1899.   u[3] = u3;
  1900.   C2length = fast_expansion_sum_zeroelim(C1length, C1, 4, u, C2);
  1901.  
  1902.   Two_Product(acxtail, bcytail, s1, s0);
  1903.   Two_Product(acytail, bcxtail, t1, t0);
  1904.   Two_Two_Diff(s1, s0, t1, t0, u3, u[2], u[1], u[0]);
  1905.   u[3] = u3;
  1906.   Dlength = fast_expansion_sum_zeroelim(C2length, C2, 4, u, D);
  1907.  
  1908.   return(D[Dlength - 1]);
  1909. }
  1910.  
  1911. double counterclockwise(
  1912. point pa,
  1913. point pb,
  1914. point pc)
  1915. {
  1916.   double detleft, detright, det;
  1917.   double detsum, errbound;
  1918.  
  1919.   counterclockcount++;
  1920.  
  1921.   detleft = (pa[0] - pc[0]) * (pb[1] - pc[1]);
  1922.   detright = (pa[1] - pc[1]) * (pb[0] - pc[0]);
  1923.   det = detleft - detright;
  1924.  
  1925.   if (detleft > 0.0) {
  1926.     if (detright <= 0.0) {
  1927.       return det;
  1928.     } else {
  1929.       detsum = detleft + detright;
  1930.     }
  1931.   } else if (detleft < 0.0) {
  1932.     if (detright >= 0.0) {
  1933.       return det;
  1934.     } else {
  1935.       detsum = -detleft - detright;
  1936.     }
  1937.   } else {
  1938.     return det;
  1939.   }
  1940.  
  1941.   errbound = ccwerrboundA * detsum;
  1942.   if ((det >= errbound) || (-det >= errbound)) {
  1943.     return det;
  1944.   }
  1945.  
  1946.   return counterclockwiseadapt(pa, pb, pc, detsum);
  1947. }
  1948.  
  1949. /*****************************************************************************/
  1950. /*                                                                           */
  1951. /*  incircle()   Return a positive value if the point pd lies inside the     */
  1952. /*               circle passing through pa, pb, and pc; a negative value if  */
  1953. /*               it lies outside; and zero if the four points are cocircular.*/
  1954. /*               The points pa, pb, and pc must be in counterclockwise       */
  1955. /*               order, or the sign of the result will be reversed.          */
  1956. /*                                                                           */
  1957. /*  Uses exact arithmetic if necessary to ensure a correct answer.  The      */
  1958. /*  result returned is the determinant of a matrix.  This determinant is     */
  1959. /*  computed adaptively, in the sense that exact arithmetic is used only to  */
  1960. /*  the degree it is needed to ensure that the returned value has the        */
  1961. /*  correct sign.  Hence, this function is usually quite fast, but will run  */
  1962. /*  more slowly when the input points are cocircular or nearly so.           */
  1963. /*                                                                           */
  1964. /*  See my Robust Predicates paper for details.                              */
  1965. /*                                                                           */
  1966. /*****************************************************************************/
  1967.  
  1968. double incircleadapt(
  1969. point pa,
  1970. point pb,
  1971. point pc,
  1972. point pd,
  1973. double permanent)
  1974. {
  1975.   INEXACT double adx, bdx, cdx, ady, bdy, cdy;
  1976.   double det, errbound;
  1977.  
  1978.   INEXACT double bdxcdy1, cdxbdy1, cdxady1, adxcdy1, adxbdy1, bdxady1;
  1979.   double bdxcdy0, cdxbdy0, cdxady0, adxcdy0, adxbdy0, bdxady0;
  1980.   double bc[4], ca[4], ab[4];
  1981.   INEXACT double bc3, ca3, ab3;
  1982.   double axbc[8], axxbc[16], aybc[8], ayybc[16], adet[32];
  1983.   int axbclen, axxbclen, aybclen, ayybclen, alen;
  1984.   double bxca[8], bxxca[16], byca[8], byyca[16], bdet[32];
  1985.   int bxcalen, bxxcalen, bycalen, byycalen, blen;
  1986.   double cxab[8], cxxab[16], cyab[8], cyyab[16], cdet[32];
  1987.   int cxablen, cxxablen, cyablen, cyyablen, clen;
  1988.   double abdet[64];
  1989.   int ablen;
  1990.   double fin1[1152], fin2[1152];
  1991.   double *finnow, *finother, *finswap;
  1992.   int finlength;
  1993.  
  1994.   double adxtail, bdxtail, cdxtail, adytail, bdytail, cdytail;
  1995.   INEXACT double adxadx1, adyady1, bdxbdx1, bdybdy1, cdxcdx1, cdycdy1;
  1996.   double adxadx0, adyady0, bdxbdx0, bdybdy0, cdxcdx0, cdycdy0;
  1997.   double aa[4], bb[4], cc[4];
  1998.   INEXACT double aa3, bb3, cc3;
  1999.   INEXACT double ti1, tj1;
  2000.   double ti0, tj0;
  2001.   double u[4], v[4];
  2002.   INEXACT double u3, v3;
  2003.   double temp8[8], temp16a[16], temp16b[16], temp16c[16];
  2004.   double temp32a[32], temp32b[32], temp48[48], temp64[64];
  2005.   int temp8len, temp16alen, temp16blen, temp16clen;
  2006.   int temp32alen, temp32blen, temp48len, temp64len;
  2007.   double axtbb[8], axtcc[8], aytbb[8], aytcc[8];
  2008.   int axtbblen, axtcclen, aytbblen, aytcclen;
  2009.   double bxtaa[8], bxtcc[8], bytaa[8], bytcc[8];
  2010.   int bxtaalen, bxtcclen, bytaalen, bytcclen;
  2011.   double cxtaa[8], cxtbb[8], cytaa[8], cytbb[8];
  2012.   int cxtaalen, cxtbblen, cytaalen, cytbblen;
  2013.   double axtbc[8], aytbc[8], bxtca[8], bytca[8], cxtab[8], cytab[8];
  2014.   int axtbclen, aytbclen, bxtcalen, bytcalen, cxtablen, cytablen;
  2015.   double axtbct[16], aytbct[16], bxtcat[16], bytcat[16], cxtabt[16], cytabt[16];
  2016.   int axtbctlen, aytbctlen, bxtcatlen, bytcatlen, cxtabtlen, cytabtlen;
  2017.   double axtbctt[8], aytbctt[8], bxtcatt[8];
  2018.   double bytcatt[8], cxtabtt[8], cytabtt[8];
  2019.   int axtbcttlen, aytbcttlen, bxtcattlen, bytcattlen, cxtabttlen, cytabttlen;
  2020.   double abt[8], bct[8], cat[8];
  2021.   int abtlen, bctlen, catlen;
  2022.   double abtt[4], bctt[4], catt[4];
  2023.   int abttlen, bcttlen, cattlen;
  2024.   INEXACT double abtt3, bctt3, catt3;
  2025.   double negate;
  2026.  
  2027.   INEXACT double bvirt;
  2028.   double avirt, bround, around;
  2029.   INEXACT double c;
  2030.   INEXACT double abig;
  2031.   double ahi, alo, bhi, blo;
  2032.   double err1, err2, err3;
  2033.   INEXACT double _i, _j;
  2034.   double _0;
  2035.  
  2036.   adx = (double) (pa[0] - pd[0]);
  2037.   bdx = (double) (pb[0] - pd[0]);
  2038.   cdx = (double) (pc[0] - pd[0]);
  2039.   ady = (double) (pa[1] - pd[1]);
  2040.   bdy = (double) (pb[1] - pd[1]);
  2041.   cdy = (double) (pc[1] - pd[1]);
  2042.  
  2043.   Two_Product(bdx, cdy, bdxcdy1, bdxcdy0);
  2044.   Two_Product(cdx, bdy, cdxbdy1, cdxbdy0);
  2045.   Two_Two_Diff(bdxcdy1, bdxcdy0, cdxbdy1, cdxbdy0, bc3, bc[2], bc[1], bc[0]);
  2046.   bc[3] = bc3;
  2047.   axbclen = scale_expansion_zeroelim(4, bc, adx, axbc);
  2048.   axxbclen = scale_expansion_zeroelim(axbclen, axbc, adx, axxbc);
  2049.   aybclen = scale_expansion_zeroelim(4, bc, ady, aybc);
  2050.   ayybclen = scale_expansion_zeroelim(aybclen, aybc, ady, ayybc);
  2051.   alen = fast_expansion_sum_zeroelim(axxbclen, axxbc, ayybclen, ayybc, adet);
  2052.  
  2053.   Two_Product(cdx, ady, cdxady1, cdxady0);
  2054.   Two_Product(adx, cdy, adxcdy1, adxcdy0);
  2055.   Two_Two_Diff(cdxady1, cdxady0, adxcdy1, adxcdy0, ca3, ca[2], ca[1], ca[0]);
  2056.   ca[3] = ca3;
  2057.   bxcalen = scale_expansion_zeroelim(4, ca, bdx, bxca);
  2058.   bxxcalen = scale_expansion_zeroelim(bxcalen, bxca, bdx, bxxca);
  2059.   bycalen = scale_expansion_zeroelim(4, ca, bdy, byca);
  2060.   byycalen = scale_expansion_zeroelim(bycalen, byca, bdy, byyca);
  2061.   blen = fast_expansion_sum_zeroelim(bxxcalen, bxxca, byycalen, byyca, bdet);
  2062.  
  2063.   Two_Product(adx, bdy, adxbdy1, adxbdy0);
  2064.   Two_Product(bdx, ady, bdxady1, bdxady0);
  2065.   Two_Two_Diff(adxbdy1, adxbdy0, bdxady1, bdxady0, ab3, ab[2], ab[1], ab[0]);
  2066.   ab[3] = ab3;
  2067.   cxablen = scale_expansion_zeroelim(4, ab, cdx, cxab);
  2068.   cxxablen = scale_expansion_zeroelim(cxablen, cxab, cdx, cxxab);
  2069.   cyablen = scale_expansion_zeroelim(4, ab, cdy, cyab);
  2070.   cyyablen = scale_expansion_zeroelim(cyablen, cyab, cdy, cyyab);
  2071.   clen = fast_expansion_sum_zeroelim(cxxablen, cxxab, cyyablen, cyyab, cdet);
  2072.  
  2073.   ablen = fast_expansion_sum_zeroelim(alen, adet, blen, bdet, abdet);
  2074.   finlength = fast_expansion_sum_zeroelim(ablen, abdet, clen, cdet, fin1);
  2075.  
  2076.   det = estimate(finlength, fin1);
  2077.   errbound = iccerrboundB * permanent;
  2078.   if ((det >= errbound) || (-det >= errbound)) {
  2079.     return det;
  2080.   }
  2081.  
  2082.   Two_Diff_Tail(pa[0], pd[0], adx, adxtail);
  2083.   Two_Diff_Tail(pa[1], pd[1], ady, adytail);
  2084.   Two_Diff_Tail(pb[0], pd[0], bdx, bdxtail);
  2085.   Two_Diff_Tail(pb[1], pd[1], bdy, bdytail);
  2086.   Two_Diff_Tail(pc[0], pd[0], cdx, cdxtail);
  2087.   Two_Diff_Tail(pc[1], pd[1], cdy, cdytail);
  2088.   if ((adxtail == 0.0) && (bdxtail == 0.0) && (cdxtail == 0.0)
  2089.       && (adytail == 0.0) && (bdytail == 0.0) && (cdytail == 0.0)) {
  2090.     return det;
  2091.   }
  2092.  
  2093.   errbound = iccerrboundC * permanent + resulterrbound * Absolute(det);
  2094.   det += ((adx * adx + ady * ady) * ((bdx * cdytail + cdy * bdxtail)
  2095.                                      - (bdy * cdxtail + cdx * bdytail))
  2096.           + 2.0 * (adx * adxtail + ady * adytail) * (bdx * cdy - bdy * cdx))
  2097.        + ((bdx * bdx + bdy * bdy) * ((cdx * adytail + ady * cdxtail)
  2098.                                      - (cdy * adxtail + adx * cdytail))
  2099.           + 2.0 * (bdx * bdxtail + bdy * bdytail) * (cdx * ady - cdy * adx))
  2100.        + ((cdx * cdx + cdy * cdy) * ((adx * bdytail + bdy * adxtail)
  2101.                                      - (ady * bdxtail + bdx * adytail))
  2102.           + 2.0 * (cdx * cdxtail + cdy * cdytail) * (adx * bdy - ady * bdx));
  2103.   if ((det >= errbound) || (-det >= errbound)) {
  2104.     return det;
  2105.   }
  2106.  
  2107.   finnow = fin1;
  2108.   finother = fin2;
  2109.  
  2110.   if ((bdxtail != 0.0) || (bdytail != 0.0)
  2111.       || (cdxtail != 0.0) || (cdytail != 0.0)) {
  2112.     Square(adx, adxadx1, adxadx0);
  2113.     Square(ady, adyady1, adyady0);
  2114.     Two_Two_Sum(adxadx1, adxadx0, adyady1, adyady0, aa3, aa[2], aa[1], aa[0]);
  2115.     aa[3] = aa3;
  2116.   }
  2117.   if ((cdxtail != 0.0) || (cdytail != 0.0)
  2118.       || (adxtail != 0.0) || (adytail != 0.0)) {
  2119.     Square(bdx, bdxbdx1, bdxbdx0);
  2120.     Square(bdy, bdybdy1, bdybdy0);
  2121.     Two_Two_Sum(bdxbdx1, bdxbdx0, bdybdy1, bdybdy0, bb3, bb[2], bb[1], bb[0]);
  2122.     bb[3] = bb3;
  2123.   }
  2124.   if ((adxtail != 0.0) || (adytail != 0.0)
  2125.       || (bdxtail != 0.0) || (bdytail != 0.0)) {
  2126.     Square(cdx, cdxcdx1, cdxcdx0);
  2127.     Square(cdy, cdycdy1, cdycdy0);
  2128.     Two_Two_Sum(cdxcdx1, cdxcdx0, cdycdy1, cdycdy0, cc3, cc[2], cc[1], cc[0]);
  2129.     cc[3] = cc3;
  2130.   }
  2131.  
  2132.   if (adxtail != 0.0) {
  2133.     axtbclen = scale_expansion_zeroelim(4, bc, adxtail, axtbc);
  2134.     temp16alen = scale_expansion_zeroelim(axtbclen, axtbc, 2.0 * adx,
  2135.                                           temp16a);
  2136.  
  2137.     axtcclen = scale_expansion_zeroelim(4, cc, adxtail, axtcc);
  2138.     temp16blen = scale_expansion_zeroelim(axtcclen, axtcc, bdy, temp16b);
  2139.  
  2140.     axtbblen = scale_expansion_zeroelim(4, bb, adxtail, axtbb);
  2141.     temp16clen = scale_expansion_zeroelim(axtbblen, axtbb, -cdy, temp16c);
  2142.  
  2143.     temp32alen = fast_expansion_sum_zeroelim(temp16alen, temp16a,
  2144.                                             temp16blen, temp16b, temp32a);
  2145.     temp48len = fast_expansion_sum_zeroelim(temp16clen, temp16c,
  2146.                                             temp32alen, temp32a, temp48);
  2147.     finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp48len,
  2148.                                             temp48, finother);
  2149.     finswap = finnow; finnow = finother; finother = finswap;
  2150.   }
  2151.   if (adytail != 0.0) {
  2152.     aytbclen = scale_expansion_zeroelim(4, bc, adytail, aytbc);
  2153.     temp16alen = scale_expansion_zeroelim(aytbclen, aytbc, 2.0 * ady,
  2154.                                           temp16a);
  2155.  
  2156.     aytbblen = scale_expansion_zeroelim(4, bb, adytail, aytbb);
  2157.     temp16blen = scale_expansion_zeroelim(aytbblen, aytbb, cdx, temp16b);
  2158.  
  2159.     aytcclen = scale_expansion_zeroelim(4, cc, adytail, aytcc);
  2160.     temp16clen = scale_expansion_zeroelim(aytcclen, aytcc, -bdx, temp16c);
  2161.  
  2162.     temp32alen = fast_expansion_sum_zeroelim(temp16alen, temp16a,
  2163.                                             temp16blen, temp16b, temp32a);
  2164.     temp48len = fast_expansion_sum_zeroelim(temp16clen, temp16c,
  2165.                                             temp32alen, temp32a, temp48);
  2166.     finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp48len,
  2167.                                             temp48, finother);
  2168.     finswap = finnow; finnow = finother; finother = finswap;
  2169.   }
  2170.   if (bdxtail != 0.0) {
  2171.     bxtcalen = scale_expansion_zeroelim(4, ca, bdxtail, bxtca);
  2172.     temp16alen = scale_expansion_zeroelim(bxtcalen, bxtca, 2.0 * bdx,
  2173.                                           temp16a);
  2174.  
  2175.     bxtaalen = scale_expansion_zeroelim(4, aa, bdxtail, bxtaa);
  2176.     temp16blen = scale_expansion_zeroelim(bxtaalen, bxtaa, cdy, temp16b);
  2177.  
  2178.     bxtcclen = scale_expansion_zeroelim(4, cc, bdxtail, bxtcc);
  2179.     temp16clen = scale_expansion_zeroelim(bxtcclen, bxtcc, -ady, temp16c);
  2180.  
  2181.     temp32alen = fast_expansion_sum_zeroelim(temp16alen, temp16a,
  2182.                                             temp16blen, temp16b, temp32a);
  2183.     temp48len = fast_expansion_sum_zeroelim(temp16clen, temp16c,
  2184.                                             temp32alen, temp32a, temp48);
  2185.     finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp48len,
  2186.                                             temp48, finother);
  2187.     finswap = finnow; finnow = finother; finother = finswap;
  2188.   }
  2189.   if (bdytail != 0.0) {
  2190.     bytcalen = scale_expansion_zeroelim(4, ca, bdytail, bytca);
  2191.     temp16alen = scale_expansion_zeroelim(bytcalen, bytca, 2.0 * bdy,
  2192.                                           temp16a);
  2193.  
  2194.     bytcclen = scale_expansion_zeroelim(4, cc, bdytail, bytcc);
  2195.     temp16blen = scale_expansion_zeroelim(bytcclen, bytcc, adx, temp16b);
  2196.  
  2197.     bytaalen = scale_expansion_zeroelim(4, aa, bdytail, bytaa);
  2198.     temp16clen = scale_expansion_zeroelim(bytaalen, bytaa, -cdx, temp16c);
  2199.  
  2200.     temp32alen = fast_expansion_sum_zeroelim(temp16alen, temp16a,
  2201.                                             temp16blen, temp16b, temp32a);
  2202.     temp48len = fast_expansion_sum_zeroelim(temp16clen, temp16c,
  2203.                                             temp32alen, temp32a, temp48);
  2204.     finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp48len,
  2205.                                             temp48, finother);
  2206.     finswap = finnow; finnow = finother; finother = finswap;
  2207.   }
  2208.   if (cdxtail != 0.0) {
  2209.     cxtablen = scale_expansion_zeroelim(4, ab, cdxtail, cxtab);
  2210.     temp16alen = scale_expansion_zeroelim(cxtablen, cxtab, 2.0 * cdx,
  2211.                                           temp16a);
  2212.  
  2213.     cxtbblen = scale_expansion_zeroelim(4, bb, cdxtail, cxtbb);
  2214.     temp16blen = scale_expansion_zeroelim(cxtbblen, cxtbb, ady, temp16b);
  2215.  
  2216.     cxtaalen = scale_expansion_zeroelim(4, aa, cdxtail, cxtaa);
  2217.     temp16clen = scale_expansion_zeroelim(cxtaalen, cxtaa, -bdy, temp16c);
  2218.  
  2219.     temp32alen = fast_expansion_sum_zeroelim(temp16alen, temp16a,
  2220.                                             temp16blen, temp16b, temp32a);
  2221.     temp48len = fast_expansion_sum_zeroelim(temp16clen, temp16c,
  2222.                                             temp32alen, temp32a, temp48);
  2223.     finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp48len,
  2224.                                             temp48, finother);
  2225.     finswap = finnow; finnow = finother; finother = finswap;
  2226.   }
  2227.   if (cdytail != 0.0) {
  2228.     cytablen = scale_expansion_zeroelim(4, ab, cdytail, cytab);
  2229.     temp16alen = scale_expansion_zeroelim(cytablen, cytab, 2.0 * cdy,
  2230.                                           temp16a);
  2231.  
  2232.     cytaalen = scale_expansion_zeroelim(4, aa, cdytail, cytaa);
  2233.     temp16blen = scale_expansion_zeroelim(cytaalen, cytaa, bdx, temp16b);
  2234.  
  2235.     cytbblen = scale_expansion_zeroelim(4, bb, cdytail, cytbb);
  2236.     temp16clen = scale_expansion_zeroelim(cytbblen, cytbb, -adx, temp16c);
  2237.  
  2238.     temp32alen = fast_expansion_sum_zeroelim(temp16alen, temp16a,
  2239.                                             temp16blen, temp16b, temp32a);
  2240.     temp48len = fast_expansion_sum_zeroelim(temp16clen, temp16c,
  2241.                                             temp32alen, temp32a, temp48);
  2242.     finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp48len,
  2243.                                             temp48, finother);
  2244.     finswap = finnow; finnow = finother; finother = finswap;
  2245.   }
  2246.  
  2247.   if ((adxtail != 0.0) || (adytail != 0.0)) {
  2248.     if ((bdxtail != 0.0) || (bdytail != 0.0)
  2249.         || (cdxtail != 0.0) || (cdytail != 0.0)) {
  2250.       Two_Product(bdxtail, cdy, ti1, ti0);
  2251.       Two_Product(bdx, cdytail, tj1, tj0);
  2252.       Two_Two_Sum(ti1, ti0, tj1, tj0, u3, u[2], u[1], u[0]);
  2253.       u[3] = u3;
  2254.       negate = -bdy;
  2255.       Two_Product(cdxtail, negate, ti1, ti0);
  2256.       negate = -bdytail;
  2257.       Two_Product(cdx, negate, tj1, tj0);
  2258.       Two_Two_Sum(ti1, ti0, tj1, tj0, v3, v[2], v[1], v[0]);
  2259.       v[3] = v3;
  2260.       bctlen = fast_expansion_sum_zeroelim(4, u, 4, v, bct);
  2261.  
  2262.       Two_Product(bdxtail, cdytail, ti1, ti0);
  2263.       Two_Product(cdxtail, bdytail, tj1, tj0);
  2264.       Two_Two_Diff(ti1, ti0, tj1, tj0, bctt3, bctt[2], bctt[1], bctt[0]);
  2265.       bctt[3] = bctt3;
  2266.       bcttlen = 4;
  2267.     } else {
  2268.       bct[0] = 0.0;
  2269.       bctlen = 1;
  2270.       bctt[0] = 0.0;
  2271.       bcttlen = 1;
  2272.     }
  2273.  
  2274.     if (adxtail != 0.0) {
  2275.       temp16alen = scale_expansion_zeroelim(axtbclen, axtbc, adxtail, temp16a);
  2276.       axtbctlen = scale_expansion_zeroelim(bctlen, bct, adxtail, axtbct);
  2277.       temp32alen = scale_expansion_zeroelim(axtbctlen, axtbct, 2.0 * adx,
  2278.                                             temp32a);
  2279.       temp48len = fast_expansion_sum_zeroelim(temp16alen, temp16a,
  2280.                                               temp32alen, temp32a, temp48);
  2281.       finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp48len,
  2282.                                               temp48, finother);
  2283.       finswap = finnow; finnow = finother; finother = finswap;
  2284.       if (bdytail != 0.0) {
  2285.         temp8len = scale_expansion_zeroelim(4, cc, adxtail, temp8);
  2286.         temp16alen = scale_expansion_zeroelim(temp8len, temp8, bdytail,
  2287.                                               temp16a);
  2288.         finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp16alen,
  2289.                                                 temp16a, finother);
  2290.         finswap = finnow; finnow = finother; finother = finswap;
  2291.       }
  2292.       if (cdytail != 0.0) {
  2293.         temp8len = scale_expansion_zeroelim(4, bb, -adxtail, temp8);
  2294.         temp16alen = scale_expansion_zeroelim(temp8len, temp8, cdytail,
  2295.                                               temp16a);
  2296.         finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp16alen,
  2297.                                                 temp16a, finother);
  2298.         finswap = finnow; finnow = finother; finother = finswap;
  2299.       }
  2300.  
  2301.       temp32alen = scale_expansion_zeroelim(axtbctlen, axtbct, adxtail,
  2302.                                             temp32a);
  2303.       axtbcttlen = scale_expansion_zeroelim(bcttlen, bctt, adxtail, axtbctt);
  2304.       temp16alen = scale_expansion_zeroelim(axtbcttlen, axtbctt, 2.0 * adx,
  2305.                                             temp16a);
  2306.       temp16blen = scale_expansion_zeroelim(axtbcttlen, axtbctt, adxtail,
  2307.                                             temp16b);
  2308.       temp32blen = fast_expansion_sum_zeroelim(temp16alen, temp16a,
  2309.                                               temp16blen, temp16b, temp32b);
  2310.       temp64len = fast_expansion_sum_zeroelim(temp32alen, temp32a,
  2311.                                               temp32blen, temp32b, temp64);
  2312.       finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp64len,
  2313.                                               temp64, finother);
  2314.       finswap = finnow; finnow = finother; finother = finswap;
  2315.     }
  2316.     if (adytail != 0.0) {
  2317.       temp16alen = scale_expansion_zeroelim(aytbclen, aytbc, adytail, temp16a);
  2318.       aytbctlen = scale_expansion_zeroelim(bctlen, bct, adytail, aytbct);
  2319.       temp32alen = scale_expansion_zeroelim(aytbctlen, aytbct, 2.0 * ady,
  2320.                                             temp32a);
  2321.       temp48len = fast_expansion_sum_zeroelim(temp16alen, temp16a,
  2322.                                               temp32alen, temp32a, temp48);
  2323.       finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp48len,
  2324.                                               temp48, finother);
  2325.       finswap = finnow; finnow = finother; finother = finswap;
  2326.  
  2327.  
  2328.       temp32alen = scale_expansion_zeroelim(aytbctlen, aytbct, adytail,
  2329.                                             temp32a);
  2330.       aytbcttlen = scale_expansion_zeroelim(bcttlen, bctt, adytail, aytbctt);
  2331.       temp16alen = scale_expansion_zeroelim(aytbcttlen, aytbctt, 2.0 * ady,
  2332.                                             temp16a);
  2333.       temp16blen = scale_expansion_zeroelim(aytbcttlen, aytbctt, adytail,
  2334.                                             temp16b);
  2335.       temp32blen = fast_expansion_sum_zeroelim(temp16alen, temp16a,
  2336.                                               temp16blen, temp16b, temp32b);
  2337.       temp64len = fast_expansion_sum_zeroelim(temp32alen, temp32a,
  2338.                                               temp32blen, temp32b, temp64);
  2339.       finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp64len,
  2340.                                               temp64, finother);
  2341.       finswap = finnow; finnow = finother; finother = finswap;
  2342.     }
  2343.   }
  2344.   if ((bdxtail != 0.0) || (bdytail != 0.0)) {
  2345.     if ((cdxtail != 0.0) || (cdytail != 0.0)
  2346.         || (adxtail != 0.0) || (adytail != 0.0)) {
  2347.       Two_Product(cdxtail, ady, ti1, ti0);
  2348.       Two_Product(cdx, adytail, tj1, tj0);
  2349.       Two_Two_Sum(ti1, ti0, tj1, tj0, u3, u[2], u[1], u[0]);
  2350.       u[3] = u3;
  2351.       negate = -cdy;
  2352.       Two_Product(adxtail, negate, ti1, ti0);
  2353.       negate = -cdytail;
  2354.       Two_Product(adx, negate, tj1, tj0);
  2355.       Two_Two_Sum(ti1, ti0, tj1, tj0, v3, v[2], v[1], v[0]);
  2356.       v[3] = v3;
  2357.       catlen = fast_expansion_sum_zeroelim(4, u, 4, v, cat);
  2358.  
  2359.       Two_Product(cdxtail, adytail, ti1, ti0);
  2360.       Two_Product(adxtail, cdytail, tj1, tj0);
  2361.       Two_Two_Diff(ti1, ti0, tj1, tj0, catt3, catt[2], catt[1], catt[0]);
  2362.       catt[3] = catt3;
  2363.       cattlen = 4;
  2364.     } else {
  2365.       cat[0] = 0.0;
  2366.       catlen = 1;
  2367.       catt[0] = 0.0;
  2368.       cattlen = 1;
  2369.     }
  2370.  
  2371.     if (bdxtail != 0.0) {
  2372.       temp16alen = scale_expansion_zeroelim(bxtcalen, bxtca, bdxtail, temp16a);
  2373.       bxtcatlen = scale_expansion_zeroelim(catlen, cat, bdxtail, bxtcat);
  2374.       temp32alen = scale_expansion_zeroelim(bxtcatlen, bxtcat, 2.0 * bdx,
  2375.                                             temp32a);
  2376.       temp48len = fast_expansion_sum_zeroelim(temp16alen, temp16a,
  2377.                                               temp32alen, temp32a, temp48);
  2378.       finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp48len,
  2379.                                               temp48, finother);
  2380.       finswap = finnow; finnow = finother; finother = finswap;
  2381.       if (cdytail != 0.0) {
  2382.         temp8len = scale_expansion_zeroelim(4, aa, bdxtail, temp8);
  2383.         temp16alen = scale_expansion_zeroelim(temp8len, temp8, cdytail,
  2384.                                               temp16a);
  2385.         finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp16alen,
  2386.                                                 temp16a, finother);
  2387.         finswap = finnow; finnow = finother; finother = finswap;
  2388.       }
  2389.       if (adytail != 0.0) {
  2390.         temp8len = scale_expansion_zeroelim(4, cc, -bdxtail, temp8);
  2391.         temp16alen = scale_expansion_zeroelim(temp8len, temp8, adytail,
  2392.                                               temp16a);
  2393.         finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp16alen,
  2394.                                                 temp16a, finother);
  2395.         finswap = finnow; finnow = finother; finother = finswap;
  2396.       }
  2397.  
  2398.       temp32alen = scale_expansion_zeroelim(bxtcatlen, bxtcat, bdxtail,
  2399.                                             temp32a);
  2400.       bxtcattlen = scale_expansion_zeroelim(cattlen, catt, bdxtail, bxtcatt);
  2401.       temp16alen = scale_expansion_zeroelim(bxtcattlen, bxtcatt, 2.0 * bdx,
  2402.                                             temp16a);
  2403.       temp16blen = scale_expansion_zeroelim(bxtcattlen, bxtcatt, bdxtail,
  2404.                                             temp16b);
  2405.       temp32blen = fast_expansion_sum_zeroelim(temp16alen, temp16a,
  2406.                                               temp16blen, temp16b, temp32b);
  2407.       temp64len = fast_expansion_sum_zeroelim(temp32alen, temp32a,
  2408.                                               temp32blen, temp32b, temp64);
  2409.       finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp64len,
  2410.                                               temp64, finother);
  2411.       finswap = finnow; finnow = finother; finother = finswap;
  2412.     }
  2413.     if (bdytail != 0.0) {
  2414.       temp16alen = scale_expansion_zeroelim(bytcalen, bytca, bdytail, temp16a);
  2415.       bytcatlen = scale_expansion_zeroelim(catlen, cat, bdytail, bytcat);
  2416.       temp32alen = scale_expansion_zeroelim(bytcatlen, bytcat, 2.0 * bdy,
  2417.                                             temp32a);
  2418.       temp48len = fast_expansion_sum_zeroelim(temp16alen, temp16a,
  2419.                                               temp32alen, temp32a, temp48);
  2420.       finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp48len,
  2421.                                               temp48, finother);
  2422.       finswap = finnow; finnow = finother; finother = finswap;
  2423.  
  2424.  
  2425.       temp32alen = scale_expansion_zeroelim(bytcatlen, bytcat, bdytail,
  2426.                                             temp32a);
  2427.       bytcattlen = scale_expansion_zeroelim(cattlen, catt, bdytail, bytcatt);
  2428.       temp16alen = scale_expansion_zeroelim(bytcattlen, bytcatt, 2.0 * bdy,
  2429.                                             temp16a);
  2430.       temp16blen = scale_expansion_zeroelim(bytcattlen, bytcatt, bdytail,
  2431.                                             temp16b);
  2432.       temp32blen = fast_expansion_sum_zeroelim(temp16alen, temp16a,
  2433.                                               temp16blen, temp16b, temp32b);
  2434.       temp64len = fast_expansion_sum_zeroelim(temp32alen, temp32a,
  2435.                                               temp32blen, temp32b, temp64);
  2436.       finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp64len,
  2437.                                               temp64, finother);
  2438.       finswap = finnow; finnow = finother; finother = finswap;
  2439.     }
  2440.   }
  2441.   if ((cdxtail != 0.0) || (cdytail != 0.0)) {
  2442.     if ((adxtail != 0.0) || (adytail != 0.0)
  2443.         || (bdxtail != 0.0) || (bdytail != 0.0)) {
  2444.       Two_Product(adxtail, bdy, ti1, ti0);
  2445.       Two_Product(adx, bdytail, tj1, tj0);
  2446.       Two_Two_Sum(ti1, ti0, tj1, tj0, u3, u[2], u[1], u[0]);
  2447.       u[3] = u3;
  2448.       negate = -ady;
  2449.       Two_Product(bdxtail, negate, ti1, ti0);
  2450.       negate = -adytail;
  2451.       Two_Product(bdx, negate, tj1, tj0);
  2452.       Two_Two_Sum(ti1, ti0, tj1, tj0, v3, v[2], v[1], v[0]);
  2453.       v[3] = v3;
  2454.       abtlen = fast_expansion_sum_zeroelim(4, u, 4, v, abt);
  2455.  
  2456.       Two_Product(adxtail, bdytail, ti1, ti0);
  2457.       Two_Product(bdxtail, adytail, tj1, tj0);
  2458.       Two_Two_Diff(ti1, ti0, tj1, tj0, abtt3, abtt[2], abtt[1], abtt[0]);
  2459.       abtt[3] = abtt3;
  2460.       abttlen = 4;
  2461.     } else {
  2462.       abt[0] = 0.0;
  2463.       abtlen = 1;
  2464.       abtt[0] = 0.0;
  2465.       abttlen = 1;
  2466.     }
  2467.  
  2468.     if (cdxtail != 0.0) {
  2469.       temp16alen = scale_expansion_zeroelim(cxtablen, cxtab, cdxtail, temp16a);
  2470.       cxtabtlen = scale_expansion_zeroelim(abtlen, abt, cdxtail, cxtabt);
  2471.       temp32alen = scale_expansion_zeroelim(cxtabtlen, cxtabt, 2.0 * cdx,
  2472.                                             temp32a);
  2473.       temp48len = fast_expansion_sum_zeroelim(temp16alen, temp16a,
  2474.                                               temp32alen, temp32a, temp48);
  2475.       finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp48len,
  2476.                                               temp48, finother);
  2477.       finswap = finnow; finnow = finother; finother = finswap;
  2478.       if (adytail != 0.0) {
  2479.         temp8len = scale_expansion_zeroelim(4, bb, cdxtail, temp8);
  2480.         temp16alen = scale_expansion_zeroelim(temp8len, temp8, adytail,
  2481.                                               temp16a);
  2482.         finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp16alen,
  2483.                                                 temp16a, finother);
  2484.         finswap = finnow; finnow = finother; finother = finswap;
  2485.       }
  2486.       if (bdytail != 0.0) {
  2487.         temp8len = scale_expansion_zeroelim(4, aa, -cdxtail, temp8);
  2488.         temp16alen = scale_expansion_zeroelim(temp8len, temp8, bdytail,
  2489.                                               temp16a);
  2490.         finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp16alen,
  2491.                                                 temp16a, finother);
  2492.         finswap = finnow; finnow = finother; finother = finswap;
  2493.       }
  2494.  
  2495.       temp32alen = scale_expansion_zeroelim(cxtabtlen, cxtabt, cdxtail,
  2496.                                             temp32a);
  2497.       cxtabttlen = scale_expansion_zeroelim(abttlen, abtt, cdxtail, cxtabtt);
  2498.       temp16alen = scale_expansion_zeroelim(cxtabttlen, cxtabtt, 2.0 * cdx,
  2499.                                             temp16a);
  2500.       temp16blen = scale_expansion_zeroelim(cxtabttlen, cxtabtt, cdxtail,
  2501.                                             temp16b);
  2502.       temp32blen = fast_expansion_sum_zeroelim(temp16alen, temp16a,
  2503.                                               temp16blen, temp16b, temp32b);
  2504.       temp64len = fast_expansion_sum_zeroelim(temp32alen, temp32a,
  2505.                                               temp32blen, temp32b, temp64);
  2506.       finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp64len,
  2507.                                               temp64, finother);
  2508.       finswap = finnow; finnow = finother; finother = finswap;
  2509.     }
  2510.     if (cdytail != 0.0) {
  2511.       temp16alen = scale_expansion_zeroelim(cytablen, cytab, cdytail, temp16a);
  2512.       cytabtlen = scale_expansion_zeroelim(abtlen, abt, cdytail, cytabt);
  2513.       temp32alen = scale_expansion_zeroelim(cytabtlen, cytabt, 2.0 * cdy,
  2514.                                             temp32a);
  2515.       temp48len = fast_expansion_sum_zeroelim(temp16alen, temp16a,
  2516.                                               temp32alen, temp32a, temp48);
  2517.       finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp48len,
  2518.                                               temp48, finother);
  2519.       finswap = finnow; finnow = finother; finother = finswap;
  2520.  
  2521.  
  2522.       temp32alen = scale_expansion_zeroelim(cytabtlen, cytabt, cdytail,
  2523.                                             temp32a);
  2524.       cytabttlen = scale_expansion_zeroelim(abttlen, abtt, cdytail, cytabtt);
  2525.       temp16alen = scale_expansion_zeroelim(cytabttlen, cytabtt, 2.0 * cdy,
  2526.                                             temp16a);
  2527.       temp16blen = scale_expansion_zeroelim(cytabttlen, cytabtt, cdytail,
  2528.                                             temp16b);
  2529.       temp32blen = fast_expansion_sum_zeroelim(temp16alen, temp16a,
  2530.                                               temp16blen, temp16b, temp32b);
  2531.       temp64len = fast_expansion_sum_zeroelim(temp32alen, temp32a,
  2532.                                               temp32blen, temp32b, temp64);
  2533.       finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp64len,
  2534.                                               temp64, finother);
  2535.       finswap = finnow; finnow = finother; finother = finswap;
  2536.     }
  2537.   }
  2538.  
  2539.   return finnow[finlength - 1];
  2540. }
  2541.  
  2542. double incircle(
  2543. point pa,
  2544. point pb,
  2545. point pc,
  2546. point pd)
  2547. {
  2548.   double adx, bdx, cdx, ady, bdy, cdy;
  2549.   double bdxcdy, cdxbdy, cdxady, adxcdy, adxbdy, bdxady;
  2550.   double alift, blift, clift;
  2551.   double det;
  2552.   double permanent, errbound;
  2553.  
  2554.   incirclecount++;
  2555.  
  2556.   adx = pa[0] - pd[0];
  2557.   bdx = pb[0] - pd[0];
  2558.   cdx = pc[0] - pd[0];
  2559.   ady = pa[1] - pd[1];
  2560.   bdy = pb[1] - pd[1];
  2561.   cdy = pc[1] - pd[1];
  2562.  
  2563.   bdxcdy = bdx * cdy;
  2564.   cdxbdy = cdx * bdy;
  2565.   alift = adx * adx + ady * ady;
  2566.  
  2567.   cdxady = cdx * ady;
  2568.   adxcdy = adx * cdy;
  2569.   blift = bdx * bdx + bdy * bdy;
  2570.  
  2571.   adxbdy = adx * bdy;
  2572.   bdxady = bdx * ady;
  2573.   clift = cdx * cdx + cdy * cdy;
  2574.  
  2575.   det = alift * (bdxcdy - cdxbdy)
  2576.       + blift * (cdxady - adxcdy)
  2577.       + clift * (adxbdy - bdxady);
  2578.  
  2579.   permanent = (Absolute(bdxcdy) + Absolute(cdxbdy)) * alift
  2580.             + (Absolute(cdxady) + Absolute(adxcdy)) * blift
  2581.             + (Absolute(adxbdy) + Absolute(bdxady)) * clift;
  2582.   errbound = iccerrboundA * permanent;
  2583.   if ((det > errbound) || (-det > errbound)) {
  2584.     return det;
  2585.   }
  2586.  
  2587.   return incircleadapt(pa, pb, pc, pd, permanent);
  2588. }
  2589.  
  2590. /**                                                                         **/
  2591. /**                                                                         **/
  2592. /********* Determinant evaluation routines end here                  *********/
  2593.  
  2594. /*****************************************************************************/
  2595. /*                                                                           */
  2596. /*  triangleinit()   Initialize some variables.                              */
  2597. /*                                                                           */
  2598. /*****************************************************************************/
  2599.  
  2600. void triangleinit(void)
  2601. {
  2602.   points.maxitems = triangles.maxitems = shelles.maxitems = viri.maxitems =
  2603.     badsegments.maxitems = badtriangles.maxitems = splaynodes.maxitems = 0l;
  2604.   points.itembytes = triangles.itembytes = shelles.itembytes = viri.itembytes =
  2605.     badsegments.itembytes = badtriangles.itembytes = splaynodes.itembytes = 0;
  2606.   recenttri.tri = (triangle *) NULL;    /* No triangle has been visited yet. */
  2607.   samples = 1;            /* Point location should take at least one sample. */
  2608.   checksegments = 0;      /* There are no segments in the triangulation yet. */
  2609.   incirclecount = counterclockcount = hyperbolacount = 0;
  2610.   circumcentercount = circletopcount = 0;
  2611.   randomseed = 1;
  2612.  
  2613.   exactinit();                     /* Initialize exact arithmetic constants. */
  2614. }
  2615.  
  2616. /*****************************************************************************/
  2617. /*                                                                           */
  2618. /*  randomnation()   Generate a random number between 0 and `choices' - 1.   */
  2619. /*                                                                           */
  2620. /*  This is a simple linear congruential random number generator.  Hence, it */
  2621. /*  is a bad random number generator, but good enough for most randomized    */
  2622. /*  geometric algorithms.                                                    */
  2623. /*                                                                           */
  2624. /*****************************************************************************/
  2625.  
  2626. unsigned long randomnation(unsigned int choices)
  2627. {
  2628.   randomseed = (randomseed * 1366l + 150889l) % 714025l;
  2629.   return randomseed / (714025l / choices + 1);
  2630. }
  2631.  
  2632.  
  2633. /********* Point location routines begin here                        *********/
  2634. /**                                                                         **/
  2635. /**                                                                         **/
  2636.  
  2637. /*****************************************************************************/
  2638. /*                                                                           */
  2639. /*  makepointmap()   Construct a mapping from points to triangles to improve  */
  2640. /*                  the speed of point location for segment insertion.       */
  2641. /*                                                                           */
  2642. /*  Traverses all the triangles, and provides each corner of each triangle   */
  2643. /*  with a pointer to that triangle.  Of course, pointers will be            */
  2644. /*  overwritten by other pointers because (almost) each point is a corner    */
  2645. /*  of several triangles, but in the end every point will point to some      */
  2646. /*  triangle that contains it.                                               */
  2647. /*                                                                           */
  2648. /*****************************************************************************/
  2649.  
  2650. void makepointmap(void)
  2651. {
  2652.   struct triedge triangleloop;
  2653.   point triorg;
  2654.  
  2655.   traversalinit(&triangles);
  2656.   triangleloop.tri = triangletraverse();
  2657.   while (triangleloop.tri != (triangle *) NULL) {
  2658.     /* Check all three points of the triangle. */
  2659.     for (triangleloop.orient = 0; triangleloop.orient < 3;
  2660.          triangleloop.orient++) {
  2661.       org(triangleloop, triorg);
  2662.       setpoint2tri(triorg, encode(triangleloop));
  2663.     }
  2664.     triangleloop.tri = triangletraverse();
  2665.   }
  2666. }
  2667.  
  2668. /*****************************************************************************/
  2669. /*                                                                           */
  2670. /*  preciselocate()   Find a triangle or edge containing a given point.      */
  2671. /*                                                                           */
  2672. /*  Begins its search from `searchtri'.  It is important that `searchtri'    */
  2673. /*  be a handle with the property that `searchpoint' is strictly to the left */
  2674. /*  of the edge denoted by `searchtri', or is collinear with that edge and   */
  2675. /*  does not intersect that edge.  (In particular, `searchpoint' should not  */
  2676. /*  be the origin or destination of that edge.)                              */
  2677. /*                                                                           */
  2678. /*  These conditions are imposed because preciselocate() is normally used in */
  2679. /*  one of two situations:                                                   */
  2680. /*                                                                           */
  2681. /*  (1)  To try to find the location to insert a new point.  Normally, we    */
  2682. /*       know an edge that the point is strictly to the left of.  In the     */
  2683. /*       incremental Delaunay algorithm, that edge is a bounding box edge.   */
  2684. /*       In Ruppert's Delaunay refinement algorithm for quality meshing,     */
  2685. /*       that edge is the shortest edge of the triangle whose circumcenter   */
  2686. /*       is being inserted.                                                  */
  2687. /*                                                                           */
  2688. /*  (2)  To try to find an existing point.  In this case, any edge on the    */
  2689. /*       convex hull is a good starting edge.  The possibility that the      */
  2690. /*       vertex one seeks is an endpoint of the starting edge must be        */
  2691. /*       screened out before preciselocate() is called.                      */
  2692. /*                                                                           */
  2693. /*  On completion, `searchtri' is a triangle that contains `searchpoint'.    */
  2694. /*                                                                           */
  2695. /*  This implementation differs from that given by Guibas and Stolfi.  It    */
  2696. /*  walks from triangle to triangle, crossing an edge only if `searchpoint'  */
  2697. /*  is on the other side of the line containing that edge.  After entering   */
  2698. /*  a triangle, there are two edges by which one can leave that triangle.    */
  2699. /*  If both edges are valid (`searchpoint' is on the other side of both      */
  2700. /*  edges), one of the two is chosen by drawing a line perpendicular to      */
  2701. /*  the entry edge (whose endpoints are `forg' and `fdest') passing through  */
  2702. /*  `fapex'.  Depending on which side of this perpendicular `searchpoint'    */
  2703. /*  falls on, an exit edge is chosen.                                        */
  2704. /*                                                                           */
  2705. /*  This implementation is empirically faster than the Guibas and Stolfi     */
  2706. /*  point location routine (which I originally used), which tends to spiral  */
  2707. /*  in toward its target.                                                    */
  2708. /*                                                                           */
  2709. /*  Returns ONVERTEX if the point lies on an existing vertex.  `searchtri'   */
  2710. /*  is a handle whose origin is the existing vertex.                         */
  2711. /*                                                                           */
  2712. /*  Returns ONEDGE if the point lies on a mesh edge.  `searchtri' is a       */
  2713. /*  handle whose primary edge is the edge on which the point lies.           */
  2714. /*                                                                           */
  2715. /*  Returns INTRIANGLE if the point lies strictly within a triangle.         */
  2716. /*  `searchtri' is a handle on the triangle that contains the point.         */
  2717. /*                                                                           */
  2718. /*  Returns OUTSIDE if the point lies outside the mesh.  `searchtri' is a    */
  2719. /*  handle whose primary edge the point is to the right of.  This might      */
  2720. /*  occur when the circumcenter of a triangle falls just slightly outside    */
  2721. /*  the mesh due to floating-point roundoff error.  It also occurs when      */
  2722. /*  seeking a hole or region point that a foolish user has placed outside    */
  2723. /*  the mesh.                                                                */
  2724. /*                                                                           */
  2725. /*  WARNING:  This routine is designed for convex triangulations, and will   */
  2726. /*  not generally work after the holes and concavities have been carved.     */
  2727. /*  However, it can still be used to find the circumcenter of a triangle, as */
  2728. /*  long as the search is begun from the triangle in question.               */
  2729. /*                                                                           */
  2730. /*****************************************************************************/
  2731.  
  2732. enum locateresult preciselocate(point searchpoint,
  2733. struct triedge *searchtri)
  2734. {
  2735.   struct triedge backtracktri;
  2736.   point forg, fdest, fapex;
  2737.   point swappoint;
  2738.   double orgorient, destorient;
  2739.   int moveleft;
  2740.   triangle ptr;                         /* Temporary variable used by sym(). */
  2741.  
  2742.   /* Where are we? */
  2743.   org(*searchtri, forg);
  2744.   dest(*searchtri, fdest);
  2745.   apex(*searchtri, fapex);
  2746.   while (1) {
  2747.     /* Check whether the apex is the point we seek. */
  2748.     if ((fapex[0] == searchpoint[0]) && (fapex[1] == searchpoint[1])) {
  2749.       lprevself(*searchtri);
  2750.       return ONVERTEX;
  2751.     }
  2752.     /* Does the point lie on the other side of the line defined by the */
  2753.     /*   triangle edge opposite the triangle's destination?            */
  2754.     destorient = counterclockwise(forg, fapex, searchpoint);
  2755.     /* Does the point lie on the other side of the line defined by the */
  2756.     /*   triangle edge opposite the triangle's origin?                 */
  2757.     orgorient = counterclockwise(fapex, fdest, searchpoint);
  2758.     if (destorient > 0.0) {
  2759.       if (orgorient > 0.0) {
  2760.         /* Move left if the inner product of (fapex - searchpoint) and  */
  2761.         /*   (fdest - forg) is positive.  This is equivalent to drawing */
  2762.         /*   a line perpendicular to the line (forg, fdest) passing     */
  2763.         /*   through `fapex', and determining which side of this line   */
  2764.         /*   `searchpoint' falls on.                                    */
  2765.         moveleft = (fapex[0] - searchpoint[0]) * (fdest[0] - forg[0]) +
  2766.                    (fapex[1] - searchpoint[1]) * (fdest[1] - forg[1]) > 0.0;
  2767.       } else {
  2768.         moveleft = 1;
  2769.       }
  2770.     } else {
  2771.       if (orgorient > 0.0) {
  2772.         moveleft = 0;
  2773.       } else {
  2774.         /* The point we seek must be on the boundary of or inside this */
  2775.         /*   triangle.                                                 */
  2776.         if (destorient == 0.0) {
  2777.           lprevself(*searchtri);
  2778.           return ONEDGE;
  2779.         }
  2780.         if (orgorient == 0.0) {
  2781.           lnextself(*searchtri);
  2782.           return ONEDGE;
  2783.         }
  2784.         return INTRIANGLE;
  2785.       }
  2786.     }
  2787.  
  2788.     /* Move to another triangle.  Leave a trace `backtracktri' in case */
  2789.     /*   floating-point roundoff or some such bogey causes us to walk  */
  2790.     /*   off a boundary of the triangulation.  We can just bounce off  */
  2791.     /*   the boundary as if it were an elastic band.                   */
  2792.     if (moveleft) {
  2793.       lprev(*searchtri, backtracktri);
  2794.       fdest = fapex;
  2795.     } else {
  2796.       lnext(*searchtri, backtracktri);
  2797.       forg = fapex;
  2798.     }
  2799.     sym(backtracktri, *searchtri);
  2800.  
  2801.     /* Check for walking off the edge. */
  2802.     if (searchtri->tri == dummytri) {
  2803.       /* Turn around. */
  2804.       triedgecopy(backtracktri, *searchtri);
  2805.       swappoint = forg;
  2806.       forg = fdest;
  2807.       fdest = swappoint;
  2808.       apex(*searchtri, fapex);
  2809.       /* Check if the point really is beyond the triangulation boundary. */
  2810.       destorient = counterclockwise(forg, fapex, searchpoint);
  2811.       orgorient = counterclockwise(fapex, fdest, searchpoint);
  2812.       if ((orgorient < 0.0) && (destorient < 0.0)) {
  2813.         return OUTSIDE;
  2814.       }
  2815.     } else {
  2816.       apex(*searchtri, fapex);
  2817.     }
  2818.   }
  2819. }
  2820.  
  2821. /*****************************************************************************/
  2822. /*                                                                           */
  2823. /*  locate()   Find a triangle or edge containing a given point.             */
  2824. /*                                                                           */
  2825. /*  Searching begins from one of:  the input `searchtri', a recently         */
  2826. /*  encountered triangle `recenttri', or from a triangle chosen from a       */
  2827. /*  random sample.  The choice is made by determining which triangle's       */
  2828. /*  origin is closest to the point we are searcing for.  Normally,           */
  2829. /*  `searchtri' should be a handle on the convex hull of the triangulation.  */
  2830. /*                                                                           */
  2831. /*  Details on the random sampling method can be found in the Mucke, Saias,  */
  2832. /*  and Zhu paper cited in the header of this code.                          */
  2833. /*                                                                           */
  2834. /*  On completion, `searchtri' is a triangle that contains `searchpoint'.    */
  2835. /*                                                                           */
  2836. /*  Returns ONVERTEX if the point lies on an existing vertex.  `searchtri'   */
  2837. /*  is a handle whose origin is the existing vertex.                         */
  2838. /*                                                                           */
  2839. /*  Returns ONEDGE if the point lies on a mesh edge.  `searchtri' is a       */
  2840. /*  handle whose primary edge is the edge on which the point lies.           */
  2841. /*                                                                           */
  2842. /*  Returns INTRIANGLE if the point lies strictly within a triangle.         */
  2843. /*  `searchtri' is a handle on the triangle that contains the point.         */
  2844. /*                                                                           */
  2845. /*  Returns OUTSIDE if the point lies outside the mesh.  `searchtri' is a    */
  2846. /*  handle whose primary edge the point is to the right of.  This might      */
  2847. /*  occur when the circumcenter of a triangle falls just slightly outside    */
  2848. /*  the mesh due to floating-point roundoff error.  It also occurs when      */
  2849. /*  seeking a hole or region point that a foolish user has placed outside    */
  2850. /*  the mesh.                                                                */
  2851. /*                                                                           */
  2852. /*  WARNING:  This routine is designed for convex triangulations, and will   */
  2853. /*  not generally work after the holes and concavities have been carved.     */
  2854. /*                                                                           */
  2855. /*****************************************************************************/
  2856.  
  2857. enum locateresult locate(point searchpoint, struct triedge *searchtri)
  2858. {
  2859.   void **sampleblock;
  2860.   triangle *firsttri;
  2861.   struct triedge sampletri;
  2862.   point torg, tdest;
  2863.   unsigned long alignptr;
  2864.   double searchdist, dist;
  2865.   double ahead;
  2866.   long sampleblocks, samplesperblock, samplenum;
  2867.   long triblocks;
  2868.   long i, j;
  2869.   triangle ptr;                         /* Temporary variable used by sym(). */
  2870.  
  2871.   /* Record the distance from the suggested starting triangle to the */
  2872.   /*   point we seek.                                                */
  2873.   org(*searchtri, torg);
  2874.   searchdist = (searchpoint[0] - torg[0]) * (searchpoint[0] - torg[0])
  2875.              + (searchpoint[1] - torg[1]) * (searchpoint[1] - torg[1]);
  2876.  
  2877.   /* If a recently encountered triangle has been recorded and has not been */
  2878.   /*   deallocated, test it as a good starting point.                      */
  2879.   if (recenttri.tri != (triangle *) NULL) {
  2880.     if (recenttri.tri[3] != (triangle) NULL) {
  2881.       org(recenttri, torg);
  2882.       if ((torg[0] == searchpoint[0]) && (torg[1] == searchpoint[1])) {
  2883.         triedgecopy(recenttri, *searchtri);
  2884.         return ONVERTEX;
  2885.       }
  2886.       dist = (searchpoint[0] - torg[0]) * (searchpoint[0] - torg[0])
  2887.            + (searchpoint[1] - torg[1]) * (searchpoint[1] - torg[1]);
  2888.       if (dist < searchdist) {
  2889.         triedgecopy(recenttri, *searchtri);
  2890.         searchdist = dist;
  2891.       }
  2892.     }
  2893.   }
  2894.  
  2895.   /* The number of random samples taken is proportional to the cube root of */
  2896.   /*   the number of triangles in the mesh.  The next bit of code assumes   */
  2897.   /*   that the number of triangles increases monotonically.                */
  2898.   while (SAMPLEFACTOR * samples * samples * samples < triangles.items) {
  2899.     samples++;
  2900.   }
  2901.   triblocks = (triangles.maxitems + TRIPERBLOCK - 1) / TRIPERBLOCK;
  2902.   samplesperblock = 1 + (samples / triblocks);
  2903.   sampleblocks = samples / samplesperblock;
  2904.   sampleblock = triangles.firstblock;
  2905.   sampletri.orient = 0;
  2906.   for (i = 0; i < sampleblocks; i++) {
  2907.     alignptr = (unsigned long) (sampleblock + 1);
  2908.     firsttri = (triangle *) (alignptr + (unsigned long) triangles.alignbytes
  2909.                           - (alignptr % (unsigned long) triangles.alignbytes));
  2910.     for (j = 0; j < samplesperblock; j++) {
  2911.       if (i == triblocks - 1) {
  2912.         samplenum = randomnation((int)
  2913.                                  (triangles.maxitems - (i * TRIPERBLOCK)));
  2914.       } else {
  2915.         samplenum = randomnation(TRIPERBLOCK);
  2916.       }
  2917.       sampletri.tri = (triangle *)
  2918.                       (firsttri + (samplenum * triangles.itemwords));
  2919.       if (sampletri.tri[3] != (triangle) NULL) {
  2920.         org(sampletri, torg);
  2921.         dist = (searchpoint[0] - torg[0]) * (searchpoint[0] - torg[0])
  2922.              + (searchpoint[1] - torg[1]) * (searchpoint[1] - torg[1]);
  2923.         if (dist < searchdist) {
  2924.           triedgecopy(sampletri, *searchtri);
  2925.           searchdist = dist;
  2926.         }
  2927.       }
  2928.     }
  2929.     sampleblock = (void **) *sampleblock;
  2930.   }
  2931.   /* Where are we? */
  2932.   org(*searchtri, torg);
  2933.   dest(*searchtri, tdest);
  2934.   /* Check the starting triangle's vertices. */
  2935.   if ((torg[0] == searchpoint[0]) && (torg[1] == searchpoint[1])) {
  2936.     return ONVERTEX;
  2937.   }
  2938.   if ((tdest[0] == searchpoint[0]) && (tdest[1] == searchpoint[1])) {
  2939.     lnextself(*searchtri);
  2940.     return ONVERTEX;
  2941.   }
  2942.   /* Orient `searchtri' to fit the preconditions of calling preciselocate(). */
  2943.   ahead = counterclockwise(torg, tdest, searchpoint);
  2944.   if (ahead < 0.0) {
  2945.     /* Turn around so that `searchpoint' is to the left of the */
  2946.     /*   edge specified by `searchtri'.                        */
  2947.     symself(*searchtri);
  2948.   } else if (ahead == 0.0) {
  2949.     /* Check if `searchpoint' is between `torg' and `tdest'. */
  2950.     if (((torg[0] < searchpoint[0]) == (searchpoint[0] < tdest[0]))
  2951.         && ((torg[1] < searchpoint[1]) == (searchpoint[1] < tdest[1]))) {
  2952.       return ONEDGE;
  2953.     }
  2954.   }
  2955.   return preciselocate(searchpoint, searchtri);
  2956. }
  2957.  
  2958. /**                                                                         **/
  2959. /**                                                                         **/
  2960. /********* Point location routines end here                          *********/
  2961.  
  2962. /********* Mesh transformation routines begin here                   *********/
  2963. /**                                                                         **/
  2964. /**                                                                         **/
  2965.  
  2966. /*****************************************************************************/
  2967. /*                                                                           */
  2968. /*  insertshelle()   Create a new shell edge and insert it between two       */
  2969. /*                   triangles.                                              */
  2970. /*                                                                           */
  2971. /*  The new shell edge is inserted at the edge described by the handle       */
  2972. /*  `tri'.  Its vertices are properly initialized.  The marker `shellemark'  */
  2973. /*  is applied to the shell edge and, if appropriate, its vertices.          */
  2974. /*                                                                           */
  2975. /*****************************************************************************/
  2976.  
  2977. void insertshelle(
  2978. struct triedge *tri,          /* Edge at which to insert the new shell edge. */
  2979. int shellemark)                            /* Marker for the new shell edge. */
  2980. {
  2981.   struct triedge oppotri;
  2982.   struct edge newshelle;
  2983.   point triorg, tridest;
  2984.   triangle ptr;                         /* Temporary variable used by sym(). */
  2985.   shelle sptr;                      /* Temporary variable used by tspivot(). */
  2986.  
  2987.   /* Mark points if possible. */
  2988.   org(*tri, triorg);
  2989.   dest(*tri, tridest);
  2990.   if (pointmark(triorg) == 0) {
  2991.     setpointmark(triorg, shellemark);
  2992.   }
  2993.   if (pointmark(tridest) == 0) {
  2994.     setpointmark(tridest, shellemark);
  2995.   }
  2996.   /* Check if there's already a shell edge here. */
  2997.   tspivot(*tri, newshelle);
  2998.   if (newshelle.sh == dummysh) {
  2999.     /* Make new shell edge and initialize its vertices. */
  3000.     makeshelle(&newshelle);
  3001.     setsorg(newshelle, tridest);
  3002.     setsdest(newshelle, triorg);
  3003.     /* Bond new shell edge to the two triangles it is sandwiched between. */
  3004.     /*   Note that the facing triangle `oppotri' might be equal to        */
  3005.     /*   `dummytri' (outer space), but the new shell edge is bonded to it */
  3006.     /*   all the same.                                                    */
  3007.     tsbond(*tri, newshelle);
  3008.     sym(*tri, oppotri);
  3009.     ssymself(newshelle);
  3010.     tsbond(oppotri, newshelle);
  3011.     setmark(newshelle, shellemark);
  3012.   } else {
  3013.     if (mark(newshelle) == 0) {
  3014.       setmark(newshelle, shellemark);
  3015.     }
  3016.   }
  3017. }
  3018.  
  3019. /*****************************************************************************/
  3020. /*                                                                           */
  3021. /*  Terminology                                                              */
  3022. /*                                                                           */
  3023. /*  A "local transformation" replaces a small set of triangles with another  */
  3024. /*  set of triangles.  This may or may not involve inserting or deleting a   */
  3025. /*  point.                                                                   */
  3026. /*                                                                           */
  3027. /*  The term "casing" is used to describe the set of triangles that are      */
  3028. /*  attached to the triangles being transformed, but are not transformed     */
  3029. /*  themselves.  Think of the casing as a fixed hollow structure inside      */
  3030. /*  which all the action happens.  A "casing" is only defined relative to    */
  3031. /*  a single transformation; each occurrence of a transformation will        */
  3032. /*  involve a different casing.                                              */
  3033. /*                                                                           */
  3034. /*  A "shell" is similar to a "casing".  The term "shell" describes the set  */
  3035. /*  of shell edges (if any) that are attached to the triangles being         */
  3036. /*  transformed.  However, I sometimes use "shell" to refer to a single      */
  3037. /*  shell edge, so don't get confused.                                       */
  3038. /*                                                                           */
  3039. /*****************************************************************************/
  3040.  
  3041. /*****************************************************************************/
  3042. /*                                                                           */
  3043. /*  flip()   Transform two triangles to two different triangles by flipping  */
  3044. /*           an edge within a quadrilateral.                                 */
  3045. /*                                                                           */
  3046. /*  Imagine the original triangles, abc and bad, oriented so that the        */
  3047. /*  shared edge ab lies in a horizontal plane, with the point b on the left  */
  3048. /*  and the point a on the right.  The point c lies below the edge, and the  */
  3049. /*  point d lies above the edge.  The `flipedge' handle holds the edge ab    */
  3050. /*  of triangle abc, and is directed left, from vertex a to vertex b.        */
  3051. /*                                                                           */
  3052. /*  The triangles abc and bad are deleted and replaced by the triangles cdb  */
  3053. /*  and dca.  The triangles that represent abc and bad are NOT deallocated;  */
  3054. /*  they are reused for dca and cdb, respectively.  Hence, any handles that  */
  3055. /*  may have held the original triangles are still valid, although not       */
  3056. /*  directed as they were before.                                            */
  3057. /*                                                                           */
  3058. /*  Upon completion of this routine, the `flipedge' handle holds the edge    */
  3059. /*  dc of triangle dca, and is directed down, from vertex d to vertex c.     */
  3060. /*  (Hence, the two triangles have rotated counterclockwise.)                */
  3061. /*                                                                           */
  3062. /*  WARNING:  This transformation is geometrically valid only if the         */
  3063. /*  quadrilateral adbc is convex.  Furthermore, this transformation is       */
  3064. /*  valid only if there is not a shell edge between the triangles abc and    */
  3065. /*  bad.  This routine does not check either of these preconditions, and     */
  3066. /*  it is the responsibility of the calling routine to ensure that they are  */
  3067. /*  met.  If they are not, the streets shall be filled with wailing and      */
  3068. /*  gnashing of teeth.                                                       */
  3069. /*                                                                           */
  3070. /*****************************************************************************/
  3071.  
  3072. void flip(
  3073. struct triedge *flipedge)                   /* Handle for the triangle abc. */
  3074. {
  3075.   struct triedge botleft, botright;
  3076.   struct triedge topleft, topright;
  3077.   struct triedge top;
  3078.   struct triedge botlcasing, botrcasing;
  3079.   struct triedge toplcasing, toprcasing;
  3080.   struct edge botlshelle, botrshelle;
  3081.   struct edge toplshelle, toprshelle;
  3082.   point leftpoint, rightpoint, botpoint;
  3083.   point farpoint;
  3084.   triangle ptr;                         /* Temporary variable used by sym(). */
  3085.   shelle sptr;                      /* Temporary variable used by tspivot(). */
  3086.  
  3087.   /* Identify the vertices of the quadrilateral. */
  3088.   org(*flipedge, rightpoint);
  3089.   dest(*flipedge, leftpoint);
  3090.   apex(*flipedge, botpoint);
  3091.   sym(*flipedge, top);
  3092.   apex(top, farpoint);
  3093.  
  3094.   /* Identify the casing of the quadrilateral. */
  3095.   lprev(top, topleft);
  3096.   sym(topleft, toplcasing);
  3097.   lnext(top, topright);
  3098.   sym(topright, toprcasing);
  3099.   lnext(*flipedge, botleft);
  3100.   sym(botleft, botlcasing);
  3101.   lprev(*flipedge, botright);
  3102.   sym(botright, botrcasing);
  3103.   /* Rotate the quadrilateral one-quarter turn counterclockwise. */
  3104.   bond(topleft, botlcasing);
  3105.   bond(botleft, botrcasing);
  3106.   bond(botright, toprcasing);
  3107.   bond(topright, toplcasing);
  3108.  
  3109.   if (checksegments) {
  3110.     /* Check for shell edges and rebond them to the quadrilateral. */
  3111.     tspivot(topleft, toplshelle);
  3112.     tspivot(botleft, botlshelle);
  3113.     tspivot(botright, botrshelle);
  3114.     tspivot(topright, toprshelle);
  3115.     if (toplshelle.sh == dummysh) {
  3116.       tsdissolve(topright);
  3117.     } else {
  3118.       tsbond(topright, toplshelle);
  3119.     }
  3120.     if (botlshelle.sh == dummysh) {
  3121.       tsdissolve(topleft);
  3122.     } else {
  3123.       tsbond(topleft, botlshelle);
  3124.     }
  3125.     if (botrshelle.sh == dummysh) {
  3126.       tsdissolve(botleft);
  3127.     } else {
  3128.       tsbond(botleft, botrshelle);
  3129.     }
  3130.     if (toprshelle.sh == dummysh) {
  3131.       tsdissolve(botright);
  3132.     } else {
  3133.       tsbond(botright, toprshelle);
  3134.     }
  3135.   }
  3136.  
  3137.   /* New point assignments for the rotated quadrilateral. */
  3138.   setorg(*flipedge, farpoint);
  3139.   setdest(*flipedge, botpoint);
  3140.   setapex(*flipedge, rightpoint);
  3141.   setorg(top, botpoint);
  3142.   setdest(top, farpoint);
  3143.   setapex(top, leftpoint);
  3144. }
  3145.  
  3146. /*****************************************************************************/
  3147. /*                                                                           */
  3148. /*  insertsite()   Insert a vertex into a Delaunay triangulation,            */
  3149. /*                 performing flips as necessary to maintain the Delaunay    */
  3150. /*                 property.                                                 */
  3151. /*                                                                           */
  3152. /*  The point `insertpoint' is located.  If `searchtri.tri' is not NULL,     */
  3153. /*  the search for the containing triangle begins from `searchtri'.  If      */
  3154. /*  `searchtri.tri' is NULL, a full point location procedure is called.      */
  3155. /*  If `insertpoint' is found inside a triangle, the triangle is split into  */
  3156. /*  three; if `insertpoint' lies on an edge, the edge is split in two,       */
  3157. /*  thereby splitting the two adjacent triangles into four.  Edge flips are  */
  3158. /*  used to restore the Delaunay property.  If `insertpoint' lies on an      */
  3159. /*  existing vertex, no action is taken, and the value DUPLICATEPOINT is     */
  3160. /*  returned.  On return, `searchtri' is set to a handle whose origin is the */
  3161. /*  existing vertex.                                                         */
  3162. /*                                                                           */
  3163. /*  Normally, the parameter `splitedge' is set to NULL, implying that no     */
  3164. /*  segment should be split.  In this case, if `insertpoint' is found to     */
  3165. /*  lie on a segment, no action is taken, and the value VIOLATINGPOINT is    */
  3166. /*  returned.  On return, `searchtri' is set to a handle whose primary edge  */
  3167. /*  is the violated segment.                                                 */
  3168. /*                                                                           */
  3169. /*  If the calling routine wishes to split a segment by inserting a point in */
  3170. /*  it, the parameter `splitedge' should be that segment.  In this case,     */
  3171. /*  `searchtri' MUST be the triangle handle reached by pivoting from that    */
  3172. /*  segment; no point location is done.                                      */
  3173. /*                                                                           */
  3174. /*  `segmentflaws' and `triflaws' are flags that indicate whether or not     */
  3175. /*  there should be checks for the creation of encroached segments or bad    */
  3176. /*  quality faces.  If a newly inserted point encroaches upon segments,      */
  3177. /*  these segments are added to the list of segments to be split if          */
  3178. /*  `segmentflaws' is set.  If bad triangles are created, these are added    */
  3179. /*  to the queue if `triflaws' is set.                                       */
  3180. /*                                                                           */
  3181. /*  If a duplicate point or violated segment does not prevent the point      */
  3182. /*  from being inserted, the return value will be ENCROACHINGPOINT if the    */
  3183. /*  point encroaches upon a segment (and checking is enabled), or            */
  3184. /*  SUCCESSFULPOINT otherwise.  In either case, `searchtri' is set to a      */
  3185. /*  handle whose origin is the newly inserted vertex.                        */
  3186. /*                                                                           */
  3187. /*  insertsite() does not use flip() for reasons of speed; some              */
  3188. /*  information can be reused from edge flip to edge flip, like the          */
  3189. /*  locations of shell edges.                                                */
  3190. /*                                                                           */
  3191. /*****************************************************************************/
  3192.  
  3193. enum insertsiteresult insertsite(
  3194. point insertpoint,
  3195. struct triedge *searchtri,
  3196. struct edge *splitedge,
  3197. int segmentflaws,
  3198. int triflaws)
  3199. {
  3200.   struct triedge horiz;
  3201.   struct triedge top;
  3202.   struct triedge botleft, botright;
  3203.   struct triedge topleft, topright;
  3204.   struct triedge newbotleft, newbotright;
  3205.   struct triedge newtopright;
  3206.   struct triedge botlcasing, botrcasing;
  3207.   struct triedge toplcasing, toprcasing;
  3208.   struct triedge testtri;
  3209.   struct edge botlshelle, botrshelle;
  3210.   struct edge toplshelle, toprshelle;
  3211.   struct edge brokenshelle;
  3212.   struct edge checkshelle;
  3213.   struct edge rightedge;
  3214.   struct edge newedge;
  3215.   struct edge *encroached;
  3216.   point first;
  3217.   point leftpoint, rightpoint, botpoint, toppoint, farpoint;
  3218.   double attrib;
  3219.   enum insertsiteresult success;
  3220.   enum locateresult intersect;
  3221.   int doflip;
  3222.   int mirrorflag;
  3223.   int i;
  3224.   triangle ptr;                         /* Temporary variable used by sym(). */
  3225.   shelle sptr;         /* Temporary variable used by spivot() and tspivot(). */
  3226.  
  3227.   if (splitedge == (struct edge *) NULL) {
  3228.     /* Find the location of the point to be inserted.  Check if a good */
  3229.     /*   starting triangle has already been provided by the caller.    */
  3230.     if (searchtri->tri == (triangle *) NULL) {
  3231.       /* Find a boundary triangle. */
  3232.       horiz.tri = dummytri;
  3233.       horiz.orient = 0;
  3234.       symself(horiz);
  3235.       /* Search for a triangle containing `insertpoint'. */
  3236.       intersect = locate(insertpoint, &horiz);
  3237.     } else {
  3238.       /* Start searching from the triangle provided by the caller. */
  3239.       triedgecopy(*searchtri, horiz);
  3240.       intersect = preciselocate(insertpoint, &horiz);
  3241.     }
  3242.   } else {
  3243.     /* The calling routine provides the edge in which the point is inserted. */
  3244.     triedgecopy(*searchtri, horiz);
  3245.     intersect = ONEDGE;
  3246.   }
  3247.   if (intersect == ONVERTEX) {
  3248.     /* There's already a vertex there.  Return in `searchtri' a triangle */
  3249.     /*   whose origin is the existing vertex.                            */
  3250.     triedgecopy(horiz, *searchtri);
  3251.     triedgecopy(horiz, recenttri);
  3252.     return DUPLICATEPOINT;
  3253.   }
  3254.   if ((intersect == ONEDGE) || (intersect == OUTSIDE)) {
  3255.     /* The vertex falls on an edge or boundary. */
  3256.     if (checksegments && (splitedge == (struct edge *) NULL)) {
  3257.       /* Check whether the vertex falls on a shell edge. */
  3258.       tspivot(horiz, brokenshelle);
  3259.       if (brokenshelle.sh != dummysh) {
  3260.         /* The vertex falls on a shell edge. */
  3261.         if (segmentflaws) {
  3262.             /* Add the shell edge to the list of encroached segments. */
  3263.             encroached = (struct edge *) poolalloc(&badsegments);
  3264.             shellecopy(brokenshelle, *encroached);
  3265.         }
  3266.         /* Return a handle whose primary edge contains the point, */
  3267.         /*   which has not been inserted.                         */
  3268.         triedgecopy(horiz, *searchtri);
  3269.         triedgecopy(horiz, recenttri);
  3270.         return VIOLATINGPOINT;
  3271.       }
  3272.     }
  3273.     /* Insert the point on an edge, dividing one triangle into two (if */
  3274.     /*   the edge lies on a boundary) or two triangles into four.      */
  3275.     lprev(horiz, botright);
  3276.     sym(botright, botrcasing);
  3277.     sym(horiz, topright);
  3278.     /* Is there a second triangle?  (Or does this edge lie on a boundary?) */
  3279.     mirrorflag = topright.tri != dummytri;
  3280.     if (mirrorflag) {
  3281.       lnextself(topright);
  3282.       sym(topright, toprcasing);
  3283.       maketriangle(&newtopright);
  3284.     } else {
  3285.       /* Splitting the boundary edge increases the number of boundary edges. */
  3286.       hullsize++;
  3287.     }
  3288.     maketriangle(&newbotright);
  3289.  
  3290.     /* Set the vertices of changed and new triangles. */
  3291.     org(horiz, rightpoint);
  3292.     dest(horiz, leftpoint);
  3293.     apex(horiz, botpoint);
  3294.     setorg(newbotright, botpoint);
  3295.     setdest(newbotright, rightpoint);
  3296.     setapex(newbotright, insertpoint);
  3297.     setorg(horiz, insertpoint);
  3298.     for (i = 0; i < eextras; i++) {
  3299.       /* Set the element attributes of a new triangle. */
  3300.       setelemattribute(newbotright, i, elemattribute(botright, i));
  3301.     }
  3302.     if (mirrorflag) {
  3303.       dest(topright, toppoint);
  3304.       setorg(newtopright, rightpoint);
  3305.       setdest(newtopright, toppoint);
  3306.       setapex(newtopright, insertpoint);
  3307.       setorg(topright, insertpoint);
  3308.       for (i = 0; i < eextras; i++) {
  3309.         /* Set the element attributes of another new triangle. */
  3310.         setelemattribute(newtopright, i, elemattribute(topright, i));
  3311.       }
  3312.     }
  3313.  
  3314.     /* There may be shell edges that need to be bonded */
  3315.     /*   to the new triangle(s).                       */
  3316.     if (checksegments) {
  3317.       tspivot(botright, botrshelle);
  3318.       if (botrshelle.sh != dummysh) {
  3319.         tsdissolve(botright);
  3320.         tsbond(newbotright, botrshelle);
  3321.       }
  3322.       if (mirrorflag) {
  3323.         tspivot(topright, toprshelle);
  3324.         if (toprshelle.sh != dummysh) {
  3325.           tsdissolve(topright);
  3326.           tsbond(newtopright, toprshelle);
  3327.         }
  3328.       }
  3329.     }
  3330.  
  3331.     /* Bond the new triangle(s) to the surrounding triangles. */
  3332.     bond(newbotright, botrcasing);
  3333.     lprevself(newbotright);
  3334.     bond(newbotright, botright);
  3335.     lprevself(newbotright);
  3336.     if (mirrorflag) {
  3337.       bond(newtopright, toprcasing);
  3338.       lnextself(newtopright);
  3339.       bond(newtopright, topright);
  3340.       lnextself(newtopright);
  3341.       bond(newtopright, newbotright);
  3342.     }
  3343.  
  3344.     if (splitedge != (struct edge *) NULL) {
  3345.       /* Split the shell edge into two. */
  3346.       setsdest(*splitedge, insertpoint);
  3347.       ssymself(*splitedge);
  3348.       spivot(*splitedge, rightedge);
  3349.       insertshelle(&newbotright, mark(*splitedge));
  3350.       tspivot(newbotright, newedge);
  3351.       sbond(*splitedge, newedge);
  3352.       ssymself(newedge);
  3353.       sbond(newedge, rightedge);
  3354.       ssymself(*splitedge);
  3355.     }
  3356.  
  3357.     /* Position `horiz' on the first edge to check for */
  3358.     /*   the Delaunay property.                        */
  3359.     lnextself(horiz);
  3360.   } else {
  3361.     /* Insert the point in a triangle, splitting it into three. */
  3362.     lnext(horiz, botleft);
  3363.     lprev(horiz, botright);
  3364.     sym(botleft, botlcasing);
  3365.     sym(botright, botrcasing);
  3366.     maketriangle(&newbotleft);
  3367.     maketriangle(&newbotright);
  3368.  
  3369.     /* Set the vertices of changed and new triangles. */
  3370.     org(horiz, rightpoint);
  3371.     dest(horiz, leftpoint);
  3372.     apex(horiz, botpoint);
  3373.     setorg(newbotleft, leftpoint);
  3374.     setdest(newbotleft, botpoint);
  3375.     setapex(newbotleft, insertpoint);
  3376.     setorg(newbotright, botpoint);
  3377.     setdest(newbotright, rightpoint);
  3378.     setapex(newbotright, insertpoint);
  3379.     setapex(horiz, insertpoint);
  3380.     for (i = 0; i < eextras; i++) {
  3381.       /* Set the element attributes of the new triangles. */
  3382.       attrib = elemattribute(horiz, i);
  3383.       setelemattribute(newbotleft, i, attrib);
  3384.       setelemattribute(newbotright, i, attrib);
  3385.     }
  3386.  
  3387.     /* There may be shell edges that need to be bonded */
  3388.     /*   to the new triangles.                         */
  3389.     if (checksegments) {
  3390.       tspivot(botleft, botlshelle);
  3391.       if (botlshelle.sh != dummysh) {
  3392.         tsdissolve(botleft);
  3393.         tsbond(newbotleft, botlshelle);
  3394.       }
  3395.       tspivot(botright, botrshelle);
  3396.       if (botrshelle.sh != dummysh) {
  3397.         tsdissolve(botright);
  3398.         tsbond(newbotright, botrshelle);
  3399.       }
  3400.     }
  3401.  
  3402.     /* Bond the new triangles to the surrounding triangles. */
  3403.     bond(newbotleft, botlcasing);
  3404.     bond(newbotright, botrcasing);
  3405.     lnextself(newbotleft);
  3406.     lprevself(newbotright);
  3407.     bond(newbotleft, newbotright);
  3408.     lnextself(newbotleft);
  3409.     bond(botleft, newbotleft);
  3410.     lprevself(newbotright);
  3411.     bond(botright, newbotright);
  3412.  
  3413.   }
  3414.  
  3415.   /* The insertion is successful by default, unless an encroached */
  3416.   /*   edge is found.                                             */
  3417.   success = SUCCESSFULPOINT;
  3418.   /* Circle around the newly inserted vertex, checking each edge opposite */
  3419.   /*   it for the Delaunay property.  Non-Delaunay edges are flipped.     */
  3420.   /*   `horiz' is always the edge being checked.  `first' marks where to  */
  3421.   /*   stop circling.                                                     */
  3422.   org(horiz, first);
  3423.   rightpoint = first;
  3424.   dest(horiz, leftpoint);
  3425.   /* Circle until finished. */
  3426.   while (1) {
  3427.     /* By default, the edge will be flipped. */
  3428.     doflip = 1;
  3429.     if (checksegments) {
  3430.       /* Check for a segment, which cannot be flipped. */
  3431.       tspivot(horiz, checkshelle);
  3432.       if (checkshelle.sh != dummysh) {
  3433.         /* The edge is a segment and cannot be flipped. */
  3434.         doflip = 0;
  3435.       }
  3436.     }
  3437.     if (doflip) {
  3438.       /* Check if the edge is a boundary edge. */
  3439.       sym(horiz, top);
  3440.       if (top.tri == dummytri) {
  3441.         /* The edge is a boundary edge and cannot be flipped. */
  3442.         doflip = 0;
  3443.       } else {
  3444.         /* Find the point on the other side of the edge. */
  3445.         apex(top, farpoint);
  3446.         /* In the incremental Delaunay triangulation algorithm, any of    */
  3447.         /*   `leftpoint', `rightpoint', and `farpoint' could be vertices  */
  3448.         /*   of the triangular bounding box.  These vertices must be      */
  3449.         /*   treated as if they are infinitely distant, even though their */
  3450.         /*   "coordinates" are not.                                       */
  3451.         if ((leftpoint == infpoint1) || (leftpoint == infpoint2)
  3452.                    || (leftpoint == infpoint3)) {
  3453.           /* `leftpoint' is infinitely distant.  Check the convexity of */
  3454.           /*   the boundary of the triangulation.  'farpoint' might be  */
  3455.           /*   infinite as well, but trust me, this same condition      */
  3456.           /*   should be applied.                                       */
  3457.           doflip = counterclockwise(insertpoint, rightpoint, farpoint) > 0.0;
  3458.         } else if ((rightpoint == infpoint1) || (rightpoint == infpoint2)
  3459.                    || (rightpoint == infpoint3)) {
  3460.           /* `rightpoint' is infinitely distant.  Check the convexity of */
  3461.           /*   the boundary of the triangulation.  'farpoint' might be  */
  3462.           /*   infinite as well, but trust me, this same condition      */
  3463.           /*   should be applied.                                       */
  3464.           doflip = counterclockwise(farpoint, leftpoint, insertpoint) > 0.0;
  3465.         } else if ((farpoint == infpoint1) || (farpoint == infpoint2)
  3466.             || (farpoint == infpoint3)) {
  3467.           /* `farpoint' is infinitely distant and cannot be inside */
  3468.           /*   the circumcircle of the triangle `horiz'.           */
  3469.           doflip = 0;
  3470.         } else {
  3471.           /* Test whether the edge is locally Delaunay. */
  3472.           doflip = incircle(leftpoint, insertpoint, rightpoint, farpoint)
  3473.                    > 0.0;
  3474.         }
  3475.         if (doflip) {
  3476.           /* We made it!  Flip the edge `horiz' by rotating its containing */
  3477.           /*   quadrilateral (the two triangles adjacent to `horiz').      */
  3478.           /* Identify the casing of the quadrilateral. */
  3479.           lprev(top, topleft);
  3480.           sym(topleft, toplcasing);
  3481.           lnext(top, topright);
  3482.           sym(topright, toprcasing);
  3483.           lnext(horiz, botleft);
  3484.           sym(botleft, botlcasing);
  3485.           lprev(horiz, botright);
  3486.           sym(botright, botrcasing);
  3487.           /* Rotate the quadrilateral one-quarter turn counterclockwise. */
  3488.           bond(topleft, botlcasing);
  3489.           bond(botleft, botrcasing);
  3490.           bond(botright, toprcasing);
  3491.           bond(topright, toplcasing);
  3492.           if (checksegments) {
  3493.             /* Check for shell edges and rebond them to the quadrilateral. */
  3494.             tspivot(topleft, toplshelle);
  3495.             tspivot(botleft, botlshelle);
  3496.             tspivot(botright, botrshelle);
  3497.             tspivot(topright, toprshelle);
  3498.             if (toplshelle.sh == dummysh) {
  3499.               tsdissolve(topright);
  3500.             } else {
  3501.               tsbond(topright, toplshelle);
  3502.             }
  3503.             if (botlshelle.sh == dummysh) {
  3504.               tsdissolve(topleft);
  3505.             } else {
  3506.               tsbond(topleft, botlshelle);
  3507.             }
  3508.             if (botrshelle.sh == dummysh) {
  3509.               tsdissolve(botleft);
  3510.             } else {
  3511.               tsbond(botleft, botrshelle);
  3512.             }
  3513.             if (toprshelle.sh == dummysh) {
  3514.               tsdissolve(botright);
  3515.             } else {
  3516.               tsbond(botright, toprshelle);
  3517.             }
  3518.           }
  3519.           /* New point assignments for the rotated quadrilateral. */
  3520.           setorg(horiz, farpoint);
  3521.           setdest(horiz, insertpoint);
  3522.           setapex(horiz, rightpoint);
  3523.           setorg(top, insertpoint);
  3524.           setdest(top, farpoint);
  3525.           setapex(top, leftpoint);
  3526.           for (i = 0; i < eextras; i++) {
  3527.             /* Take the average of the two triangles' attributes. */
  3528.             attrib = 0.5 * (elemattribute(top, i) + elemattribute(horiz, i));
  3529.             setelemattribute(top, i, attrib);
  3530.             setelemattribute(horiz, i, attrib);
  3531.           }
  3532.           /* On the next iterations, consider the two edges that were  */
  3533.           /*   exposed (this is, are now visible to the newly inserted */
  3534.           /*   point) by the edge flip.                                */
  3535.           lprevself(horiz);
  3536.           leftpoint = farpoint;
  3537.         }
  3538.       }
  3539.     }
  3540.     if (!doflip) {
  3541.       /* The handle `horiz' is accepted as locally Delaunay. */
  3542.       /* Look for the next edge around the newly inserted point. */
  3543.       lnextself(horiz);
  3544.       sym(horiz, testtri);
  3545.       /* Check for finishing a complete revolution about the new point, or */
  3546.       /*   falling off the edge of the triangulation.  The latter will     */
  3547.       /*   happen when a point is inserted at a boundary.                  */
  3548.       if ((leftpoint == first) || (testtri.tri == dummytri)) {
  3549.         /* We're done.  Return a triangle whose origin is the new point. */
  3550.         lnext(horiz, *searchtri);
  3551.         lnext(horiz, recenttri);
  3552.         return success;
  3553.       }
  3554.       /* Finish finding the next edge around the newly inserted point. */
  3555.       lnext(testtri, horiz);
  3556.       rightpoint = leftpoint;
  3557.       dest(horiz, leftpoint);
  3558.     }
  3559.   }
  3560. }
  3561.  
  3562. /*****************************************************************************/
  3563. /*                                                                           */
  3564. /*  triangulatepolygon()   Find the Delaunay triangulation of a polygon that */
  3565. /*                         has a certain "nice" shape.  This includes the    */
  3566. /*                         polygons that result from deletion of a point or  */
  3567. /*                         insertion of a segment.                           */
  3568. /*                                                                           */
  3569. /*  This is a conceptually difficult routine.  The starting assumption is    */
  3570. /*  that we have a polygon with n sides.  n - 1 of these sides are currently */
  3571. /*  represented as edges in the mesh.  One side, called the "base", need not */
  3572. /*  be.                                                                      */
  3573. /*                                                                           */
  3574. /*  Inside the polygon is a structure I call a "fan", consisting of n - 1    */
  3575. /*  triangles that share a common origin.  For each of these triangles, the  */
  3576. /*  edge opposite the origin is one of the sides of the polygon.  The        */
  3577. /*  primary edge of each triangle is the edge directed from the origin to    */
  3578. /*  the destination; note that this is not the same edge that is a side of   */
  3579. /*  the polygon.  `firstedge' is the primary edge of the first triangle.     */
  3580. /*  From there, the triangles follow in counterclockwise order about the     */
  3581. /*  polygon, until `lastedge', the primary edge of the last triangle.        */
  3582. /*  `firstedge' and `lastedge' are probably connected to other triangles     */
  3583. /*  beyond the extremes of the fan, but their identity is not important, as  */
  3584. /*  long as the fan remains connected to them.                               */
  3585. /*                                                                           */
  3586. /*  Imagine the polygon oriented so that its base is at the bottom.  This    */
  3587. /*  puts `firstedge' on the far right, and `lastedge' on the far left.       */
  3588. /*  The right vertex of the base is the destination of `firstedge', and the  */
  3589. /*  left vertex of the base is the apex of `lastedge'.                       */
  3590. /*                                                                           */
  3591. /*  The challenge now is to find the right sequence of edge flips to         */
  3592. /*  transform the fan into a Delaunay triangulation of the polygon.  Each    */
  3593. /*  edge flip effectively removes one triangle from the fan, committing it   */
  3594. /*  to the polygon.  The resulting polygon has one fewer edge.  If `doflip'  */
  3595. /*  is set, the final flip will be performed, resulting in a fan of one      */
  3596. /*  (useless?) triangle.  If `doflip' is not set, the final flip is not      */
  3597. /*  performed, resulting in a fan of two triangles, and an unfinished        */
  3598. /*  triangular polygon that is not yet filled out with a single triangle.    */
  3599. /*  On completion of the routine, `lastedge' is the last remaining triangle, */
  3600. /*  or the leftmost of the last two.                                         */
  3601. /*                                                                           */
  3602. /*  Although the flips are performed in the order described above, the       */
  3603. /*  decisions about what flips to perform are made in precisely the reverse  */
  3604. /*  order.  The recursive triangulatepolygon() procedure makes a decision,   */
  3605. /*  uses up to two recursive calls to triangulate the "subproblems"          */
  3606. /*  (polygons with fewer edges), and then performs an edge flip.             */
  3607. /*                                                                           */
  3608. /*  The "decision" it makes is which vertex of the polygon should be         */
  3609. /*  connected to the base.  This decision is made by testing every possible  */
  3610. /*  vertex.  Once the best vertex is found, the two edges that connect this  */
  3611. /*  vertex to the base become the bases for two smaller polygons.  These     */
  3612. /*  are triangulated recursively.  Unfortunately, this approach can take     */
  3613. /*  O(n^2) time not only in the worst case, but in many common cases.  It's  */
  3614. /*  rarely a big deal for point deletion, where n is rarely larger than ten, */
  3615. /*  but it could be a big deal for segment insertion, especially if there's  */
  3616. /*  a lot of long segments that each cut many triangles.  I ought to code    */
  3617. /*  a faster algorithm some time.                                            */
  3618. /*                                                                           */
  3619. /*  The `edgecount' parameter is the number of sides of the polygon,         */
  3620. /*  including its base.  `triflaws' is a flag that determines whether the    */
  3621. /*  new triangles should be tested for quality, and enqueued if they are     */
  3622. /*  bad.                                                                     */
  3623. /*                                                                           */
  3624. /*****************************************************************************/
  3625.  
  3626. void triangulatepolygon(
  3627. struct triedge *firstedge,
  3628. struct triedge *lastedge,
  3629. int edgecount,
  3630. int doflip,
  3631. int triflaws)
  3632. {
  3633.   struct triedge testtri;
  3634.   struct triedge besttri;
  3635.   struct triedge tempedge;
  3636.   point leftbasepoint, rightbasepoint;
  3637.   point testpoint;
  3638.   point bestpoint;
  3639.   int bestnumber;
  3640.   int i;
  3641.   triangle ptr;   /* Temporary variable used by sym(), onext(), and oprev(). */
  3642.  
  3643.   /* Identify the base vertices. */
  3644.   apex(*lastedge, leftbasepoint);
  3645.   dest(*firstedge, rightbasepoint);
  3646.   /* Find the best vertex to connect the base to. */
  3647.   onext(*firstedge, besttri);
  3648.   dest(besttri, bestpoint);
  3649.   triedgecopy(besttri, testtri);
  3650.   bestnumber = 1;
  3651.   for (i = 2; i <= edgecount - 2; i++) {
  3652.     onextself(testtri);
  3653.     dest(testtri, testpoint);
  3654.     /* Is this a better vertex? */
  3655.     if (incircle(leftbasepoint, rightbasepoint, bestpoint, testpoint) > 0.0) {
  3656.       triedgecopy(testtri, besttri);
  3657.       bestpoint = testpoint;
  3658.       bestnumber = i;
  3659.     }
  3660.   }
  3661.   if (bestnumber > 1) {
  3662.     /* Recursively triangulate the smaller polygon on the right. */
  3663.     oprev(besttri, tempedge);
  3664.     triangulatepolygon(firstedge, &tempedge, bestnumber + 1, 1, triflaws);
  3665.   }
  3666.   if (bestnumber < edgecount - 2) {
  3667.     /* Recursively triangulate the smaller polygon on the left. */
  3668.     sym(besttri, tempedge);
  3669.     triangulatepolygon(&besttri, lastedge, edgecount - bestnumber, 1,
  3670.                        triflaws);
  3671.     /* Find `besttri' again; it may have been lost to edge flips. */
  3672.     sym(tempedge, besttri);
  3673.   }
  3674.   if (doflip) {
  3675.     /* Do one final edge flip. */
  3676.     flip(&besttri);
  3677.   }
  3678.   /* Return the base triangle. */
  3679.   triedgecopy(besttri, *lastedge);
  3680. }
  3681.  
  3682.  
  3683. /**                                                                         **/
  3684. /**                                                                         **/
  3685. /********* Mesh transformation routines end here                     *********/
  3686.  
  3687. /********* Divide-and-conquer Delaunay triangulation begins here     *********/
  3688. /**                                                                         **/
  3689. /**                                                                         **/
  3690.  
  3691. /*****************************************************************************/
  3692. /*                                                                           */
  3693. /*  The divide-and-conquer bounding box                                      */
  3694. /*                                                                           */
  3695. /*  I originally implemented the divide-and-conquer and incremental Delaunay */
  3696. /*  triangulations using the edge-based data structure presented by Guibas   */
  3697. /*  and Stolfi.  Switching to a triangle-based data structure doubled the    */
  3698. /*  speed.  However, I had to think of a few extra tricks to maintain the    */
  3699. /*  elegance of the original algorithms.                                     */
  3700. /*                                                                           */
  3701. /*  The "bounding box" used by my variant of the divide-and-conquer          */
  3702. /*  algorithm uses one triangle for each edge of the convex hull of the      */
  3703. /*  triangulation.  These bounding triangles all share a common apical       */
  3704. /*  vertex, which is represented by NULL and which represents nothing.       */
  3705. /*  The bounding triangles are linked in a circular fan about this NULL      */
  3706. /*  vertex, and the edges on the convex hull of the triangulation appear     */
  3707. /*  opposite the NULL vertex.  You might find it easiest to imagine that     */
  3708. /*  the NULL vertex is a point in 3D space behind the center of the          */
  3709. /*  triangulation, and that the bounding triangles form a sort of cone.      */
  3710. /*                                                                           */
  3711. /*  This bounding box makes it easy to represent degenerate cases.  For      */
  3712. /*  instance, the triangulation of two vertices is a single edge.  This edge */
  3713. /*  is represented by two bounding box triangles, one on each "side" of the  */
  3714. /*  edge.  These triangles are also linked together in a fan about the NULL  */
  3715. /*  vertex.                                                                  */
  3716. /*                                                                           */
  3717. /*  The bounding box also makes it easy to traverse the convex hull, as the  */
  3718. /*  divide-and-conquer algorithm needs to do.                                */
  3719. /*                                                                           */
  3720. /*****************************************************************************/
  3721.  
  3722. /*****************************************************************************/
  3723. /*                                                                           */
  3724. /*  pointsort()   Sort an array of points by x-coordinate, using the         */
  3725. /*                y-coordinate as a secondary key.                           */
  3726. /*                                                                           */
  3727. /*  Uses quicksort.  Randomized O(n log n) time.  No, I did not make any of  */
  3728. /*  the usual quicksort mistakes.                                            */
  3729. /*                                                                           */
  3730. /*****************************************************************************/
  3731.  
  3732. void pointsort(
  3733. point *sortarray,
  3734. int arraysize)
  3735. {
  3736.   int left, right;
  3737.   int pivot;
  3738.   double pivotx, pivoty;
  3739.   point temp;
  3740.  
  3741.   if (arraysize == 2) {
  3742.     /* Recursive base case. */
  3743.     if ((sortarray[0][0] > sortarray[1][0]) ||
  3744.         ((sortarray[0][0] == sortarray[1][0]) &&
  3745.          (sortarray[0][1] > sortarray[1][1]))) {
  3746.       temp = sortarray[1];
  3747.       sortarray[1] = sortarray[0];
  3748.       sortarray[0] = temp;
  3749.     }
  3750.     return;
  3751.   }
  3752.   /* Choose a random pivot to split the array. */
  3753.   pivot = (int) randomnation(arraysize);
  3754.   pivotx = sortarray[pivot][0];
  3755.   pivoty = sortarray[pivot][1];
  3756.   /* Split the array. */
  3757.   left = -1;
  3758.   right = arraysize;
  3759.   while (left < right) {
  3760.     /* Search for a point whose x-coordinate is too large for the left. */
  3761.     do {
  3762.       left++;
  3763.     } while ((left <= right) && ((sortarray[left][0] < pivotx) ||
  3764.                                  ((sortarray[left][0] == pivotx) &&
  3765.                                   (sortarray[left][1] < pivoty))));
  3766.     /* Search for a point whose x-coordinate is too small for the right. */
  3767.     do {
  3768.       right--;
  3769.     } while ((left <= right) && ((sortarray[right][0] > pivotx) ||
  3770.                                  ((sortarray[right][0] == pivotx) &&
  3771.                                   (sortarray[right][1] > pivoty))));
  3772.     if (left < right) {
  3773.       /* Swap the left and right points. */
  3774.       temp = sortarray[left];
  3775.       sortarray[left] = sortarray[right];
  3776.       sortarray[right] = temp;
  3777.     }
  3778.   }
  3779.   if (left > 1) {
  3780.     /* Recursively sort the left subset. */
  3781.     pointsort(sortarray, left);
  3782.   }
  3783.   if (right < arraysize - 2) {
  3784.     /* Recursively sort the right subset. */
  3785.     pointsort(&sortarray[right + 1], arraysize - right - 1);
  3786.   }
  3787. }
  3788.  
  3789. /*****************************************************************************/
  3790. /*                                                                           */
  3791. /*  pointmedian()   An order statistic algorithm, almost.  Shuffles an array */
  3792. /*                  of points so that the first `median' points occur        */
  3793. /*                  lexicographically before the remaining points.           */
  3794. /*                                                                           */
  3795. /*  Uses the x-coordinate as the primary key if axis == 0; the y-coordinate  */
  3796. /*  if axis == 1.  Very similar to the pointsort() procedure, but runs in    */
  3797. /*  randomized linear time.                                                  */
  3798. /*                                                                           */
  3799. /*****************************************************************************/
  3800.  
  3801. void pointmedian(
  3802. point *sortarray,
  3803. int arraysize,
  3804. int median,
  3805. int axis)
  3806. {
  3807.   int left, right;
  3808.   int pivot;
  3809.   double pivot1, pivot2;
  3810.   point temp;
  3811.  
  3812.   if (arraysize == 2) {
  3813.     /* Recursive base case. */
  3814.     if ((sortarray[0][axis] > sortarray[1][axis]) ||
  3815.         ((sortarray[0][axis] == sortarray[1][axis]) &&
  3816.          (sortarray[0][1 - axis] > sortarray[1][1 - axis]))) {
  3817.       temp = sortarray[1];
  3818.       sortarray[1] = sortarray[0];
  3819.       sortarray[0] = temp;
  3820.     }
  3821.     return;
  3822.   }
  3823.   /* Choose a random pivot to split the array. */
  3824.   pivot = (int) randomnation(arraysize);
  3825.   pivot1 = sortarray[pivot][axis];
  3826.   pivot2 = sortarray[pivot][1 - axis];
  3827.   /* Split the array. */
  3828.   left = -1;
  3829.   right = arraysize;
  3830.   while (left < right) {
  3831.     /* Search for a point whose x-coordinate is too large for the left. */
  3832.     do {
  3833.       left++;
  3834.     } while ((left <= right) && ((sortarray[left][axis] < pivot1) ||
  3835.                                  ((sortarray[left][axis] == pivot1) &&
  3836.                                   (sortarray[left][1 - axis] < pivot2))));
  3837.     /* Search for a point whose x-coordinate is too small for the right. */
  3838.     do {
  3839.       right--;
  3840.     } while ((left <= right) && ((sortarray[right][axis] > pivot1) ||
  3841.                                  ((sortarray[right][axis] == pivot1) &&
  3842.                                   (sortarray[right][1 - axis] > pivot2))));
  3843.     if (left < right) {
  3844.       /* Swap the left and right points. */
  3845.       temp = sortarray[left];
  3846.       sortarray[left] = sortarray[right];
  3847.       sortarray[right] = temp;
  3848.     }
  3849.   }
  3850.   /* Unlike in pointsort(), at most one of the following */
  3851.   /*   conditionals is true.                             */
  3852.   if (left > median) {
  3853.     /* Recursively shuffle the left subset. */
  3854.     pointmedian(sortarray, left, median, axis);
  3855.   }
  3856.   if (right < median - 1) {
  3857.     /* Recursively shuffle the right subset. */
  3858.     pointmedian(&sortarray[right + 1], arraysize - right - 1,
  3859.                 median - right - 1, axis);
  3860.   }
  3861. }
  3862.  
  3863. /*****************************************************************************/
  3864. /*                                                                           */
  3865. /*  alternateaxes()   Sorts the points as appropriate for the divide-and-    */
  3866. /*                    conquer algorithm with alternating cuts.               */
  3867. /*                                                                           */
  3868. /*  Partitions by x-coordinate if axis == 0; by y-coordinate if axis == 1.   */
  3869. /*  For the base case, subsets containing only two or three points are       */
  3870. /*  always sorted by x-coordinate.                                           */
  3871. /*                                                                           */
  3872. /*****************************************************************************/
  3873.  
  3874. void alternateaxes(point *sortarray,
  3875. int arraysize,
  3876. int axis)
  3877. {
  3878.   int divider;
  3879.  
  3880.   divider = arraysize >> 1;
  3881.   if (arraysize <= 3) {
  3882.     /* Recursive base case:  subsets of two or three points will be      */
  3883.     /*   handled specially, and should always be sorted by x-coordinate. */
  3884.     axis = 0;
  3885.   }
  3886.   /* Partition with a horizontal or vertical cut. */
  3887.   pointmedian(sortarray, arraysize, divider, axis);
  3888.   /* Recursively partition the subsets with a cross cut. */
  3889.   if (arraysize - divider >= 2) {
  3890.     if (divider >= 2) {
  3891.       alternateaxes(sortarray, divider, 1 - axis);
  3892.     }
  3893.     alternateaxes(&sortarray[divider], arraysize - divider, 1 - axis);
  3894.   }
  3895. }
  3896.  
  3897. /*****************************************************************************/
  3898. /*                                                                           */
  3899. /*  mergehulls()   Merge two adjacent Delaunay triangulations into a         */
  3900. /*                 single Delaunay triangulation.                            */
  3901. /*                                                                           */
  3902. /*  This is similar to the algorithm given by Guibas and Stolfi, but uses    */
  3903. /*  a triangle-based, rather than edge-based, data structure.                */
  3904. /*                                                                           */
  3905. /*  The algorithm walks up the gap between the two triangulations, knitting  */
  3906. /*  them together.  As they are merged, some of their bounding triangles     */
  3907. /*  are converted into real triangles of the triangulation.  The procedure   */
  3908. /*  pulls each hull's bounding triangles apart, then knits them together     */
  3909. /*  like the teeth of two gears.  The Delaunay property determines, at each  */
  3910. /*  step, whether the next "tooth" is a bounding triangle of the left hull   */
  3911. /*  or the right.  When a bounding triangle becomes real, its apex is        */
  3912. /*  changed from NULL to a real point.                                       */
  3913. /*                                                                           */
  3914. /*  Only two new triangles need to be allocated.  These become new bounding  */
  3915. /*  triangles at the top and bottom of the seam.  They are used to connect   */
  3916. /*  the remaining bounding triangles (those that have not been converted     */
  3917. /*  into real triangles) into a single fan.                                  */
  3918. /*                                                                           */
  3919. /*  On entry, `farleft' and `innerleft' are bounding triangles of the left   */
  3920. /*  triangulation.  The origin of `farleft' is the leftmost vertex, and      */
  3921. /*  the destination of `innerleft' is the rightmost vertex of the            */
  3922. /*  triangulation.  Similarly, `innerright' and `farright' are bounding      */
  3923. /*  triangles of the right triangulation.  The origin of `innerright' and    */
  3924. /*  destination of `farright' are the leftmost and rightmost vertices.       */
  3925. /*                                                                           */
  3926. /*  On completion, the origin of `farleft' is the leftmost vertex of the     */
  3927. /*  merged triangulation, and the destination of `farright' is the rightmost */
  3928. /*  vertex.                                                                  */
  3929. /*                                                                           */
  3930. /*****************************************************************************/
  3931.  
  3932. void mergehulls(
  3933. struct triedge *farleft,
  3934. struct triedge *innerleft,
  3935. struct triedge *innerright,
  3936. struct triedge *farright,
  3937. int axis)
  3938. {
  3939.   struct triedge leftcand, rightcand;
  3940.   struct triedge baseedge;
  3941.   struct triedge nextedge;
  3942.   struct triedge sidecasing, topcasing, outercasing;
  3943.   struct triedge checkedge;
  3944.   point innerleftdest;
  3945.   point innerrightorg;
  3946.   point innerleftapex, innerrightapex;
  3947.   point farleftpt, farrightpt;
  3948.   point farleftapex, farrightapex;
  3949.   point lowerleft, lowerright;
  3950.   point upperleft, upperright;
  3951.   point nextapex;
  3952.   point checkvertex;
  3953.   int changemade;
  3954.   int badedge;
  3955.   int leftfinished, rightfinished;
  3956.   triangle ptr;                         /* Temporary variable used by sym(). */
  3957.  
  3958.   dest(*innerleft, innerleftdest);
  3959.   apex(*innerleft, innerleftapex);
  3960.   org(*innerright, innerrightorg);
  3961.   apex(*innerright, innerrightapex);
  3962.   /* Special treatment for horizontal cuts. */
  3963.   if (axis == 1) {
  3964.     org(*farleft, farleftpt);
  3965.     apex(*farleft, farleftapex);
  3966.     dest(*farright, farrightpt);
  3967.     apex(*farright, farrightapex);
  3968.     /* The pointers to the extremal points are shifted to point to the */
  3969.     /*   topmost and bottommost point of each hull, rather than the    */
  3970.     /*   leftmost and rightmost points.                                */
  3971.     while (farleftapex[1] < farleftpt[1]) {
  3972.       lnextself(*farleft);
  3973.       symself(*farleft);
  3974.       farleftpt = farleftapex;
  3975.       apex(*farleft, farleftapex);
  3976.     }
  3977.     sym(*innerleft, checkedge);
  3978.     apex(checkedge, checkvertex);
  3979.     while (checkvertex[1] > innerleftdest[1]) {
  3980.       lnext(checkedge, *innerleft);
  3981.       innerleftapex = innerleftdest;
  3982.       innerleftdest = checkvertex;
  3983.       sym(*innerleft, checkedge);
  3984.       apex(checkedge, checkvertex);
  3985.     }
  3986.     while (innerrightapex[1] < innerrightorg[1]) {
  3987.       lnextself(*innerright);
  3988.       symself(*innerright);
  3989.       innerrightorg = innerrightapex;
  3990.       apex(*innerright, innerrightapex);
  3991.     }
  3992.     sym(*farright, checkedge);
  3993.     apex(checkedge, checkvertex);
  3994.     while (checkvertex[1] > farrightpt[1]) {
  3995.       lnext(checkedge, *farright);
  3996.       farrightapex = farrightpt;
  3997.       farrightpt = checkvertex;
  3998.       sym(*farright, checkedge);
  3999.       apex(checkedge, checkvertex);
  4000.     }
  4001.   }
  4002.   /* Find a line tangent to and below both hulls. */
  4003.   do {
  4004.     changemade = 0;
  4005.     /* Make innerleftdest the "bottommost" point of the left hull. */
  4006.     if (counterclockwise(innerleftdest, innerleftapex, innerrightorg) > 0.0) {
  4007.       lprevself(*innerleft);
  4008.       symself(*innerleft);
  4009.       innerleftdest = innerleftapex;
  4010.       apex(*innerleft, innerleftapex);
  4011.       changemade = 1;
  4012.     }
  4013.     /* Make innerrightorg the "bottommost" point of the right hull. */
  4014.     if (counterclockwise(innerrightapex, innerrightorg, innerleftdest) > 0.0) {
  4015.       lnextself(*innerright);
  4016.       symself(*innerright);
  4017.       innerrightorg = innerrightapex;
  4018.       apex(*innerright, innerrightapex);
  4019.       changemade = 1;
  4020.     }
  4021.   } while (changemade);
  4022.   /* Find the two candidates to be the next "gear tooth". */
  4023.   sym(*innerleft, leftcand);
  4024.   sym(*innerright, rightcand);
  4025.   /* Create the bottom new bounding triangle. */
  4026.   maketriangle(&baseedge);
  4027.   /* Connect it to the bounding boxes of the left and right triangulations. */
  4028.   bond(baseedge, *innerleft);
  4029.   lnextself(baseedge);
  4030.   bond(baseedge, *innerright);
  4031.   lnextself(baseedge);
  4032.   setorg(baseedge, innerrightorg);
  4033.   setdest(baseedge, innerleftdest);
  4034.   /* Apex is intentionally left NULL. */
  4035.   /* Fix the extreme triangles if necessary. */
  4036.   org(*farleft, farleftpt);
  4037.   if (innerleftdest == farleftpt) {
  4038.     lnext(baseedge, *farleft);
  4039.   }
  4040.   dest(*farright, farrightpt);
  4041.   if (innerrightorg == farrightpt) {
  4042.     lprev(baseedge, *farright);
  4043.   }
  4044.   /* The vertices of the current knitting edge. */
  4045.   lowerleft = innerleftdest;
  4046.   lowerright = innerrightorg;
  4047.   /* The candidate vertices for knitting. */
  4048.   apex(leftcand, upperleft);
  4049.   apex(rightcand, upperright);
  4050.   /* Walk up the gap between the two triangulations, knitting them together. */
  4051.   while (1) {
  4052.     /* Have we reached the top?  (This isn't quite the right question,       */
  4053.     /*   because even though the left triangulation might seem finished now, */
  4054.     /*   moving up on the right triangulation might reveal a new point of    */
  4055.     /*   the left triangulation.  And vice-versa.)                           */
  4056.     leftfinished = counterclockwise(upperleft, lowerleft, lowerright) <= 0.0;
  4057.     rightfinished = counterclockwise(upperright, lowerleft, lowerright) <= 0.0;
  4058.     if (leftfinished && rightfinished) {
  4059.       /* Create the top new bounding triangle. */
  4060.       maketriangle(&nextedge);
  4061.       setorg(nextedge, lowerleft);
  4062.       setdest(nextedge, lowerright);
  4063.       /* Apex is intentionally left NULL. */
  4064.       /* Connect it to the bounding boxes of the two triangulations. */
  4065.       bond(nextedge, baseedge);
  4066.       lnextself(nextedge);
  4067.       bond(nextedge, rightcand);
  4068.       lnextself(nextedge);
  4069.       bond(nextedge, leftcand);
  4070.       /* Special treatment for horizontal cuts. */
  4071.       if (axis == 1) {
  4072.         org(*farleft, farleftpt);
  4073.         apex(*farleft, farleftapex);
  4074.         dest(*farright, farrightpt);
  4075.         apex(*farright, farrightapex);
  4076.         sym(*farleft, checkedge);
  4077.         apex(checkedge, checkvertex);
  4078.         /* The pointers to the extremal points are restored to the leftmost */
  4079.         /*   and rightmost points (rather than topmost and bottommost).     */
  4080.         while (checkvertex[0] < farleftpt[0]) {
  4081.           lprev(checkedge, *farleft);
  4082.           farleftapex = farleftpt;
  4083.           farleftpt = checkvertex;
  4084.           sym(*farleft, checkedge);
  4085.           apex(checkedge, checkvertex);
  4086.         }
  4087.         while (farrightapex[0] > farrightpt[0]) {
  4088.           lprevself(*farright);
  4089.           symself(*farright);
  4090.           farrightpt = farrightapex;
  4091.           apex(*farright, farrightapex);
  4092.         }
  4093.       }
  4094.       return;
  4095.     }
  4096.     /* Consider eliminating edges from the left triangulation. */
  4097.     if (!leftfinished) {
  4098.       /* What vertex would be exposed if an edge were deleted? */
  4099.       lprev(leftcand, nextedge);
  4100.       symself(nextedge);
  4101.       apex(nextedge, nextapex);
  4102.       /* If nextapex is NULL, then no vertex would be exposed; the */
  4103.       /*   triangulation would have been eaten right through.      */
  4104.       if (nextapex != (point) NULL) {
  4105.         /* Check whether the edge is Delaunay. */
  4106.         badedge = incircle(lowerleft, lowerright, upperleft, nextapex) > 0.0;
  4107.         while (badedge) {
  4108.           /* Eliminate the edge with an edge flip.  As a result, the    */
  4109.           /*   left triangulation will have one more boundary triangle. */
  4110.           lnextself(nextedge);
  4111.           sym(nextedge, topcasing);
  4112.           lnextself(nextedge);
  4113.           sym(nextedge, sidecasing);
  4114.           bond(nextedge, topcasing);
  4115.           bond(leftcand, sidecasing);
  4116.           lnextself(leftcand);
  4117.           sym(leftcand, outercasing);
  4118.           lprevself(nextedge);
  4119.           bond(nextedge, outercasing);
  4120.           /* Correct the vertices to reflect the edge flip. */
  4121.           setorg(leftcand, lowerleft);
  4122.           setdest(leftcand, NULL);
  4123.           setapex(leftcand, nextapex);
  4124.           setorg(nextedge, NULL);
  4125.           setdest(nextedge, upperleft);
  4126.           setapex(nextedge, nextapex);
  4127.           /* Consider the newly exposed vertex. */
  4128.           upperleft = nextapex;
  4129.           /* What vertex would be exposed if another edge were deleted? */
  4130.           triedgecopy(sidecasing, nextedge);
  4131.           apex(nextedge, nextapex);
  4132.           if (nextapex != (point) NULL) {
  4133.             /* Check whether the edge is Delaunay. */
  4134.             badedge = incircle(lowerleft, lowerright, upperleft, nextapex)
  4135.                       > 0.0;
  4136.           } else {
  4137.             /* Avoid eating right through the triangulation. */
  4138.             badedge = 0;
  4139.           }
  4140.         }
  4141.       }
  4142.     }
  4143.     /* Consider eliminating edges from the right triangulation. */
  4144.     if (!rightfinished) {
  4145.       /* What vertex would be exposed if an edge were deleted? */
  4146.       lnext(rightcand, nextedge);
  4147.       symself(nextedge);
  4148.       apex(nextedge, nextapex);
  4149.       /* If nextapex is NULL, then no vertex would be exposed; the */
  4150.       /*   triangulation would have been eaten right through.      */
  4151.       if (nextapex != (point) NULL) {
  4152.         /* Check whether the edge is Delaunay. */
  4153.         badedge = incircle(lowerleft, lowerright, upperright, nextapex) > 0.0;
  4154.         while (badedge) {
  4155.           /* Eliminate the edge with an edge flip.  As a result, the     */
  4156.           /*   right triangulation will have one more boundary triangle. */
  4157.           lprevself(nextedge);
  4158.           sym(nextedge, topcasing);
  4159.           lprevself(nextedge);
  4160.           sym(nextedge, sidecasing);
  4161.           bond(nextedge, topcasing);
  4162.           bond(rightcand, sidecasing);
  4163.           lprevself(rightcand);
  4164.           sym(rightcand, outercasing);
  4165.           lnextself(nextedge);
  4166.           bond(nextedge, outercasing);
  4167.           /* Correct the vertices to reflect the edge flip. */
  4168.           setorg(rightcand, NULL);
  4169.           setdest(rightcand, lowerright);
  4170.           setapex(rightcand, nextapex);
  4171.           setorg(nextedge, upperright);
  4172.           setdest(nextedge, NULL);
  4173.           setapex(nextedge, nextapex);
  4174.           /* Consider the newly exposed vertex. */
  4175.           upperright = nextapex;
  4176.           /* What vertex would be exposed if another edge were deleted? */
  4177.           triedgecopy(sidecasing, nextedge);
  4178.           apex(nextedge, nextapex);
  4179.           if (nextapex != (point) NULL) {
  4180.             /* Check whether the edge is Delaunay. */
  4181.             badedge = incircle(lowerleft, lowerright, upperright, nextapex)
  4182.                       > 0.0;
  4183.           } else {
  4184.             /* Avoid eating right through the triangulation. */
  4185.             badedge = 0;
  4186.           }
  4187.         }
  4188.       }
  4189.     }
  4190.     if (leftfinished || (!rightfinished &&
  4191.            (incircle(upperleft, lowerleft, lowerright, upperright) > 0.0))) {
  4192.       /* Knit the triangulations, adding an edge from `lowerleft' */
  4193.       /*   to `upperright'.                                       */
  4194.       bond(baseedge, rightcand);
  4195.       lprev(rightcand, baseedge);
  4196.       setdest(baseedge, lowerleft);
  4197.       lowerright = upperright;
  4198.       sym(baseedge, rightcand);
  4199.       apex(rightcand, upperright);
  4200.     } else {
  4201.       /* Knit the triangulations, adding an edge from `upperleft' */
  4202.       /*   to `lowerright'.                                       */
  4203.       bond(baseedge, leftcand);
  4204.       lnext(leftcand, baseedge);
  4205.       setorg(baseedge, lowerright);
  4206.       lowerleft = upperleft;
  4207.       sym(baseedge, leftcand);
  4208.       apex(leftcand, upperleft);
  4209.     }
  4210.   }
  4211. }
  4212.  
  4213. /*****************************************************************************/
  4214. /*                                                                           */
  4215. /*  divconqrecurse()   Recursively form a Delaunay triangulation by the      */
  4216. /*                     divide-and-conquer method.                            */
  4217. /*                                                                           */
  4218. /*  Recursively breaks down the problem into smaller pieces, which are       */
  4219. /*  knitted together by mergehulls().  The base cases (problems of two or    */
  4220. /*  three points) are handled specially here.                                */
  4221. /*                                                                           */
  4222. /*  On completion, `farleft' and `farright' are bounding triangles such that */
  4223. /*  the origin of `farleft' is the leftmost vertex (breaking ties by         */
  4224. /*  choosing the highest leftmost vertex), and the destination of            */
  4225. /*  `farright' is the rightmost vertex (breaking ties by choosing the        */
  4226. /*  lowest rightmost vertex).                                                */
  4227. /*                                                                           */
  4228. /*****************************************************************************/
  4229.  
  4230. void divconqrecurse(
  4231. point *sortarray,
  4232. int vertices,
  4233. int axis,
  4234. struct triedge *farleft,
  4235. struct triedge *farright)
  4236. {
  4237.   struct triedge midtri, tri1, tri2, tri3;
  4238.   struct triedge innerleft, innerright;
  4239.   double area;
  4240.   int divider;
  4241.  
  4242.   if (vertices == 2) {
  4243.     /* The triangulation of two vertices is an edge.  An edge is */
  4244.     /*   represented by two bounding triangles.                  */
  4245.     maketriangle(farleft);
  4246.     setorg(*farleft, sortarray[0]);
  4247.     setdest(*farleft, sortarray[1]);
  4248.     /* The apex is intentionally left NULL. */
  4249.     maketriangle(farright);
  4250.     setorg(*farright, sortarray[1]);
  4251.     setdest(*farright, sortarray[0]);
  4252.     /* The apex is intentionally left NULL. */
  4253.     bond(*farleft, *farright);
  4254.     lprevself(*farleft);
  4255.     lnextself(*farright);
  4256.     bond(*farleft, *farright);
  4257.     lprevself(*farleft);
  4258.     lnextself(*farright);
  4259.     bond(*farleft, *farright);
  4260.     /* Ensure that the origin of `farleft' is sortarray[0]. */
  4261.     lprev(*farright, *farleft);
  4262.     return;
  4263.   } else if (vertices == 3) {
  4264.     /* The triangulation of three vertices is either a triangle (with */
  4265.     /*   three bounding triangles) or two edges (with four bounding   */
  4266.     /*   triangles).  In either case, four triangles are created.     */
  4267.     maketriangle(&midtri);
  4268.     maketriangle(&tri1);
  4269.     maketriangle(&tri2);
  4270.     maketriangle(&tri3);
  4271.     area = counterclockwise(sortarray[0], sortarray[1], sortarray[2]);
  4272.     if (area == 0.0) {
  4273.       /* Three collinear points; the triangulation is two edges. */
  4274.       setorg(midtri, sortarray[0]);
  4275.       setdest(midtri, sortarray[1]);
  4276.       setorg(tri1, sortarray[1]);
  4277.       setdest(tri1, sortarray[0]);
  4278.       setorg(tri2, sortarray[2]);
  4279.       setdest(tri2, sortarray[1]);
  4280.       setorg(tri3, sortarray[1]);
  4281.       setdest(tri3, sortarray[2]);
  4282.       /* All apices are intentionally left NULL. */
  4283.       bond(midtri, tri1);
  4284.       bond(tri2, tri3);
  4285.       lnextself(midtri);
  4286.       lprevself(tri1);
  4287.       lnextself(tri2);
  4288.       lprevself(tri3);
  4289.       bond(midtri, tri3);
  4290.       bond(tri1, tri2);
  4291.       lnextself(midtri);
  4292.       lprevself(tri1);
  4293.       lnextself(tri2);
  4294.       lprevself(tri3);
  4295.       bond(midtri, tri1);
  4296.       bond(tri2, tri3);
  4297.       /* Ensure that the origin of `farleft' is sortarray[0]. */
  4298.       triedgecopy(tri1, *farleft);
  4299.       /* Ensure that the destination of `farright' is sortarray[2]. */
  4300.       triedgecopy(tri2, *farright);
  4301.     } else {
  4302.       /* The three points are not collinear; the triangulation is one */
  4303.       /*   triangle, namely `midtri'.                                 */
  4304.       setorg(midtri, sortarray[0]);
  4305.       setdest(tri1, sortarray[0]);
  4306.       setorg(tri3, sortarray[0]);
  4307.       /* Apices of tri1, tri2, and tri3 are left NULL. */
  4308.       if (area > 0.0) {
  4309.         /* The vertices are in counterclockwise order. */
  4310.         setdest(midtri, sortarray[1]);
  4311.         setorg(tri1, sortarray[1]);
  4312.         setdest(tri2, sortarray[1]);
  4313.         setapex(midtri, sortarray[2]);
  4314.         setorg(tri2, sortarray[2]);
  4315.         setdest(tri3, sortarray[2]);
  4316.       } else {
  4317.         /* The vertices are in clockwise order. */
  4318.         setdest(midtri, sortarray[2]);
  4319.         setorg(tri1, sortarray[2]);
  4320.         setdest(tri2, sortarray[2]);
  4321.         setapex(midtri, sortarray[1]);
  4322.         setorg(tri2, sortarray[1]);
  4323.         setdest(tri3, sortarray[1]);
  4324.       }
  4325.       /* The topology does not depend on how the vertices are ordered. */
  4326.       bond(midtri, tri1);
  4327.       lnextself(midtri);
  4328.       bond(midtri, tri2);
  4329.       lnextself(midtri);
  4330.       bond(midtri, tri3);
  4331.       lprevself(tri1);
  4332.       lnextself(tri2);
  4333.       bond(tri1, tri2);
  4334.       lprevself(tri1);
  4335.       lprevself(tri3);
  4336.       bond(tri1, tri3);
  4337.       lnextself(tri2);
  4338.       lprevself(tri3);
  4339.       bond(tri2, tri3);
  4340.       /* Ensure that the origin of `farleft' is sortarray[0]. */
  4341.       triedgecopy(tri1, *farleft);
  4342.       /* Ensure that the destination of `farright' is sortarray[2]. */
  4343.       if (area > 0.0) {
  4344.         triedgecopy(tri2, *farright);
  4345.       } else {
  4346.         lnext(*farleft, *farright);
  4347.       }
  4348.     }
  4349.     return;
  4350.   } else {
  4351.     /* Split the vertices in half. */
  4352.     divider = vertices >> 1;
  4353.     /* Recursively triangulate each half. */
  4354.     divconqrecurse(sortarray, divider, 1 - axis, farleft, &innerleft);
  4355.     divconqrecurse(&sortarray[divider], vertices - divider, 1 - axis,
  4356.                    &innerright, farright);
  4357.     /* Merge the two triangulations into one. */
  4358.     mergehulls(farleft, &innerleft, &innerright, farright, axis);
  4359.   }
  4360. }
  4361.  
  4362. long removeghosts(struct triedge *startghost)
  4363. {
  4364.   struct triedge searchedge;
  4365.   struct triedge dissolveedge;
  4366.   struct triedge deadtri;
  4367.   long hullsize;
  4368.   triangle ptr;                         /* Temporary variable used by sym(). */
  4369.  
  4370.   /* Find an edge on the convex hull to start point location from. */
  4371.   lprev(*startghost, searchedge);
  4372.   symself(searchedge);
  4373.   dummytri[0] = encode(searchedge);
  4374.   /* Remove the bounding box and count the convex hull edges. */
  4375.   triedgecopy(*startghost, dissolveedge);
  4376.   hullsize = 0;
  4377.   do {
  4378.     hullsize++;
  4379.     lnext(dissolveedge, deadtri);
  4380.     lprevself(dissolveedge);
  4381.     symself(dissolveedge);
  4382.     /* Remove a bounding triangle from a convex hull triangle. */
  4383.     dissolve(dissolveedge);
  4384.     /* Find the next bounding triangle. */
  4385.     sym(deadtri, dissolveedge);
  4386.     /* Delete the bounding triangle. */
  4387.     triangledealloc(deadtri.tri);
  4388.   } while (!triedgeequal(dissolveedge, *startghost));
  4389.   return hullsize;
  4390. }
  4391.  
  4392. /*****************************************************************************/
  4393. /*                                                                           */
  4394. /*  divconqdelaunay()   Form a Delaunay triangulation by the divide-and-     */
  4395. /*                      conquer method.                                      */
  4396. /*                                                                           */
  4397. /*  Sorts the points, calls a recursive procedure to triangulate them, and   */
  4398. /*  removes the bounding box, setting boundary markers as appropriate.       */
  4399. /*                                                                           */
  4400. /*****************************************************************************/
  4401.  
  4402. long divconqdelaunay(void)
  4403. {
  4404.   point *sortarray;
  4405.   struct triedge hullleft, hullright;
  4406.   int divider;
  4407.   int i, j;
  4408.  
  4409.   /* Allocate an array of pointers to points for sorting. */
  4410.   sortarray = (point *) malloc(inpoints * sizeof(point));
  4411.   if (sortarray == (point *) NULL) {
  4412.     vTrace("Error:  Out of memory.");
  4413.     exit(1);
  4414.   }
  4415.   traversalinit(&points);
  4416.   for (i = 0; i < inpoints; i++) {
  4417.     sortarray[i] = pointtraverse();
  4418.   }
  4419.   /* Sort the points. */
  4420.   pointsort(sortarray, inpoints);
  4421.   /* Discard duplicate points, which can really mess up the algorithm. */
  4422.   i = 0;
  4423.   for (j = 1; j < inpoints; j++) {
  4424.     if ((sortarray[i][0] == sortarray[j][0])
  4425.         && (sortarray[i][1] == sortarray[j][1])) {
  4426. /*  Commented out - would eliminate point from output .node file, but causes
  4427.     a failure if some segment has this point as an endpoint.
  4428.       setpointmark(sortarray[j], DEADPOINT);
  4429. */
  4430.     } else {
  4431.       i++;
  4432.       sortarray[i] = sortarray[j];
  4433.     }
  4434.   }
  4435.   i++;
  4436.     /* Re-sort the array of points to accommodate alternating cuts. */
  4437.     divider = i >> 1;
  4438.     if (i - divider >= 2) {
  4439.       if (divider >= 2) {
  4440.         alternateaxes(sortarray, divider, 1);
  4441.       }
  4442.       alternateaxes(&sortarray[divider], i - divider, 1);
  4443.     }
  4444.   /* Form the Delaunay triangulation. */
  4445.   divconqrecurse(sortarray, i, 0, &hullleft, &hullright);
  4446.   free(sortarray);
  4447.  
  4448.   return removeghosts(&hullleft);
  4449. }
  4450.  
  4451. /**                                                                         **/
  4452. /**                                                                         **/
  4453. /********* Divide-and-conquer Delaunay triangulation ends here       *********/
  4454.  
  4455. /********* General mesh construction routines begin here             *********/
  4456. /**                                                                         **/
  4457. /**                                                                         **/
  4458.  
  4459. /*****************************************************************************/
  4460. /*                                                                           */
  4461. /*  delaunay()   Form a Delaunay triangulation.                              */
  4462. /*                                                                           */
  4463. /*****************************************************************************/
  4464.  
  4465. long delaunay(void)
  4466. {
  4467.   eextras = 0;
  4468.   initializetrisegpools();
  4469.  
  4470.   return divconqdelaunay();
  4471. }
  4472.  
  4473. /**                                                                         **/
  4474. /**                                                                         **/
  4475. /********* General mesh construction routines end here               *********/
  4476.  
  4477. /********* Segment (shell edge) insertion begins here                *********/
  4478. /**                                                                         **/
  4479. /**                                                                         **/
  4480.  
  4481. /*****************************************************************************/
  4482. /*                                                                           */
  4483. /*  finddirection()   Find the first triangle on the path from one point     */
  4484. /*                    to another.                                            */
  4485. /*                                                                           */
  4486. /*  Finds the triangle that intersects a line segment drawn from the         */
  4487. /*  origin of `searchtri' to the point `endpoint', and returns the result    */
  4488. /*  in `searchtri'.  The origin of `searchtri' does not change, even though  */
  4489. /*  the triangle returned may differ from the one passed in.  This routine   */
  4490. /*  is used to find the direction to move in to get from one point to        */
  4491. /*  another.                                                                 */
  4492. /*                                                                           */
  4493. /*  The return value notes whether the destination or apex of the found      */
  4494. /*  triangle is collinear with the two points in question.                   */
  4495. /*                                                                           */
  4496. /*****************************************************************************/
  4497.  
  4498. enum finddirectionresult finddirection(
  4499. struct triedge *searchtri,
  4500. point endpoint)
  4501. {
  4502.   struct triedge checktri;
  4503.   point startpoint;
  4504.   point leftpoint, rightpoint;
  4505.   double leftccw, rightccw;
  4506.   int leftflag, rightflag;
  4507.   triangle ptr;           /* Temporary variable used by onext() and oprev(). */
  4508.  
  4509.   org(*searchtri, startpoint);
  4510.   dest(*searchtri, rightpoint);
  4511.   apex(*searchtri, leftpoint);
  4512.   /* Is `endpoint' to the left? */
  4513.   leftccw = counterclockwise(endpoint, startpoint, leftpoint);
  4514.   leftflag = leftccw > 0.0;
  4515.   /* Is `endpoint' to the right? */
  4516.   rightccw = counterclockwise(startpoint, endpoint, rightpoint);
  4517.   rightflag = rightccw > 0.0;
  4518.   if (leftflag && rightflag) {
  4519.     /* `searchtri' faces directly away from `endpoint'.  We could go */
  4520.     /*   left or right.  Ask whether it's a triangle or a boundary   */
  4521.     /*   on the left.                                                */
  4522.     onext(*searchtri, checktri);
  4523.     if (checktri.tri == dummytri) {
  4524.       leftflag = 0;
  4525.     } else {
  4526.       rightflag = 0;
  4527.     }
  4528.   }
  4529.   while (leftflag) {
  4530.     /* Turn left until satisfied. */
  4531.     onextself(*searchtri);
  4532.     if (searchtri->tri == dummytri) {
  4533.       vTrace("Internal error in finddirection():  Unable to find a");
  4534.       vTrace("  triangle leading from (%.12g, %.12g) to", startpoint[0],
  4535.              startpoint[1]);
  4536.       vTrace("  (%.12g, %.12g).", endpoint[0], endpoint[1]);
  4537.       internalerror();
  4538.     }
  4539.     apex(*searchtri, leftpoint);
  4540.     rightccw = leftccw;
  4541.     leftccw = counterclockwise(endpoint, startpoint, leftpoint);
  4542.     leftflag = leftccw > 0.0;
  4543.   }
  4544.   while (rightflag) {
  4545.     /* Turn right until satisfied. */
  4546.     oprevself(*searchtri);
  4547.     if (searchtri->tri == dummytri) {
  4548.       vTrace("Internal error in finddirection():  Unable to find a");
  4549.       vTrace("  triangle leading from (%.12g, %.12g) to", startpoint[0],
  4550.              startpoint[1]);
  4551.       vTrace("  (%.12g, %.12g).", endpoint[0], endpoint[1]);
  4552.       internalerror();
  4553.     }
  4554.     dest(*searchtri, rightpoint);
  4555.     leftccw = rightccw;
  4556.     rightccw = counterclockwise(startpoint, endpoint, rightpoint);
  4557.     rightflag = rightccw > 0.0;
  4558.   }
  4559.   if (leftccw == 0.0) {
  4560.     return LEFTCOLLINEAR;
  4561.   } else if (rightccw == 0.0) {
  4562.     return RIGHTCOLLINEAR;
  4563.   } else {
  4564.     return WITHIN;
  4565.   }
  4566. }
  4567.  
  4568. /*****************************************************************************/
  4569. /*                                                                           */
  4570. /*  segmentintersection()   Find the intersection of an existing segment     */
  4571. /*                          and a segment that is being inserted.  Insert    */
  4572. /*                          a point at the intersection, splitting an        */
  4573. /*                          existing shell edge.                             */
  4574. /*                                                                           */
  4575. /*  The segment being inserted connects the apex of splittri to endpoint2.   */
  4576. /*  splitshelle is the shell edge being split, and MUST be opposite          */
  4577. /*  splittri.  Hence, the edge being split connects the origin and           */
  4578. /*  destination of splittri.                                                 */
  4579. /*                                                                           */
  4580. /*  On completion, splittri is a handle having the newly inserted            */
  4581. /*  intersection point as its origin, and endpoint1 as its destination.      */
  4582. /*                                                                           */
  4583. /*****************************************************************************/
  4584.  
  4585. void segmentintersection(
  4586. struct triedge *splittri,
  4587. struct edge *splitshelle,
  4588. point endpoint2)
  4589. {
  4590.   point endpoint1;
  4591.   point torg, tdest;
  4592.   point leftpoint, rightpoint;
  4593.   point newpoint;
  4594.   enum insertsiteresult success;
  4595.   enum finddirectionresult collinear;
  4596.   double ex, ey;
  4597.   double tx, ty;
  4598.   double etx, ety;
  4599.   double split, denom;
  4600.   int i;
  4601.   triangle ptr;                       /* Temporary variable used by onext(). */
  4602.  
  4603.   /* Find the other three segment endpoints. */
  4604.   apex(*splittri, endpoint1);
  4605.   org(*splittri, torg);
  4606.   dest(*splittri, tdest);
  4607.   /* Segment intersection formulae; see the Antonio reference. */
  4608.   tx = tdest[0] - torg[0];
  4609.   ty = tdest[1] - torg[1];
  4610.   ex = endpoint2[0] - endpoint1[0];
  4611.   ey = endpoint2[1] - endpoint1[1];
  4612.   etx = torg[0] - endpoint2[0];
  4613.   ety = torg[1] - endpoint2[1];
  4614.   denom = ty * ex - tx * ey;
  4615.   if (denom == 0.0) {
  4616.     vTrace("Internal error in segmentintersection():");
  4617.     vTrace("  Attempt to find intersection of parallel segments.");
  4618.     internalerror();
  4619.   }
  4620.   split = (ey * etx - ex * ety) / denom;
  4621.   /* Create the new point. */
  4622.   newpoint = (point) poolalloc(&points);
  4623.   /* Interpolate its coordinate and attributes. */
  4624.   for (i = 0; i < 2 + nextras; i++) {
  4625.     newpoint[i] = torg[i] + split * (tdest[i] - torg[i]);
  4626.   }
  4627.   setpointmark(newpoint, mark(*splitshelle));
  4628.   /* Insert the intersection point.  This should always succeed. */
  4629.   success = insertsite(newpoint, splittri, splitshelle, 0, 0);
  4630.   if (success != SUCCESSFULPOINT) {
  4631.     vTrace("Internal error in segmentintersection():");
  4632.     vTrace("  Failure to split a segment.");
  4633.     internalerror();
  4634.   }
  4635.   /* Inserting the point may have caused edge flips.  We wish to rediscover */
  4636.   /*   the edge connecting endpoint1 to the new intersection point.         */
  4637.   collinear = finddirection(splittri, endpoint1);
  4638.   dest(*splittri, rightpoint);
  4639.   apex(*splittri, leftpoint);
  4640.   if ((leftpoint[0] == endpoint1[0]) && (leftpoint[1] == endpoint1[1])) {
  4641.     onextself(*splittri);
  4642.   } else if ((rightpoint[0] != endpoint1[0]) ||
  4643.              (rightpoint[1] != endpoint1[1])) {
  4644.     vTrace("Internal error in segmentintersection():");
  4645.     vTrace("  Topological inconsistency after splitting a segment.");
  4646.     internalerror();
  4647.   }
  4648.   /* `splittri' should have destination endpoint1. */
  4649. }
  4650.  
  4651. /*****************************************************************************/
  4652. /*                                                                           */
  4653. /*  scoutsegment()   Scout the first triangle on the path from one endpoint  */
  4654. /*                   to another, and check for completion (reaching the      */
  4655. /*                   second endpoint), a collinear point, and the            */
  4656. /*                   intersection of two segments.                           */
  4657. /*                                                                           */
  4658. /*  Returns one if the entire segment is successfully inserted, and zero if  */
  4659. /*  the job must be finished by conformingedge() or constrainededge().       */
  4660. /*                                                                           */
  4661. /*  If the first triangle on the path has the second endpoint as its         */
  4662. /*  destination or apex, a shell edge is inserted and the job is done.       */
  4663. /*                                                                           */
  4664. /*  If the first triangle on the path has a destination or apex that lies on */
  4665. /*  the segment, a shell edge is inserted connecting the first endpoint to   */
  4666. /*  the collinear point, and the search is continued from the collinear      */
  4667. /*  point.                                                                   */
  4668. /*                                                                           */
  4669. /*  If the first triangle on the path has a shell edge opposite its origin,  */
  4670. /*  then there is a segment that intersects the segment being inserted.      */
  4671. /*  Their intersection point is inserted, splitting the shell edge.          */
  4672. /*                                                                           */
  4673. /*  Otherwise, return zero.                                                  */
  4674. /*                                                                           */
  4675. /*****************************************************************************/
  4676.  
  4677. int scoutsegment(
  4678. struct triedge *searchtri,
  4679. point endpoint2,
  4680. int newmark)
  4681. {
  4682.   struct triedge crosstri;
  4683.   struct edge crossedge;
  4684.   point leftpoint, rightpoint;
  4685.   point endpoint1;
  4686.   enum finddirectionresult collinear;
  4687.   shelle sptr;                      /* Temporary variable used by tspivot(). */
  4688.  
  4689.   collinear = finddirection(searchtri, endpoint2);
  4690.   dest(*searchtri, rightpoint);
  4691.   apex(*searchtri, leftpoint);
  4692.   if (((leftpoint[0] == endpoint2[0]) && (leftpoint[1] == endpoint2[1])) ||
  4693.       ((rightpoint[0] == endpoint2[0]) && (rightpoint[1] == endpoint2[1]))) {
  4694.     /* The segment is already an edge in the mesh. */
  4695.     if ((leftpoint[0] == endpoint2[0]) && (leftpoint[1] == endpoint2[1])) {
  4696.       lprevself(*searchtri);
  4697.     }
  4698.     /* Insert a shell edge, if there isn't already one there. */
  4699.     insertshelle(searchtri, newmark);
  4700.     return 1;
  4701.   } else if (collinear == LEFTCOLLINEAR) {
  4702.     /* We've collided with a point between the segment's endpoints. */
  4703.     /* Make the collinear point be the triangle's origin. */
  4704.     lprevself(*searchtri);
  4705.     insertshelle(searchtri, newmark);
  4706.     /* Insert the remainder of the segment. */
  4707.     return scoutsegment(searchtri, endpoint2, newmark);
  4708.   } else if (collinear == RIGHTCOLLINEAR) {
  4709.     /* We've collided with a point between the segment's endpoints. */
  4710.     insertshelle(searchtri, newmark);
  4711.     /* Make the collinear point be the triangle's origin. */
  4712.     lnextself(*searchtri);
  4713.     /* Insert the remainder of the segment. */
  4714.     return scoutsegment(searchtri, endpoint2, newmark);
  4715.   } else {
  4716.     lnext(*searchtri, crosstri);
  4717.     tspivot(crosstri, crossedge);
  4718.     /* Check for a crossing segment. */
  4719.     if (crossedge.sh == dummysh) {
  4720.       return 0;
  4721.     } else {
  4722.       org(*searchtri, endpoint1);
  4723.       /* Insert a point at the intersection. */
  4724.       segmentintersection(&crosstri, &crossedge, endpoint2);
  4725.       triedgecopy(crosstri, *searchtri);
  4726.       insertshelle(searchtri, newmark);
  4727.       /* Insert the remainder of the segment. */
  4728.       return scoutsegment(searchtri, endpoint2, newmark);
  4729.     }
  4730.   }
  4731. }
  4732.  
  4733. /*****************************************************************************/
  4734. /*                                                                           */
  4735. /*  delaunayfixup()   Enforce the Delaunay condition at an edge, fanning out */
  4736. /*                    recursively from an existing point.  Pay special       */
  4737. /*                    attention to stacking inverted triangles.              */
  4738. /*                                                                           */
  4739. /*  This is a support routine for inserting segments into a constrained      */
  4740. /*  Delaunay triangulation.                                                  */
  4741. /*                                                                           */
  4742. /*  The origin of fixuptri is treated as if it has just been inserted, and   */
  4743. /*  the local Delaunay condition needs to be enforced.  It is only enforced  */
  4744. /*  in one sector, however, that being the angular range defined by          */
  4745. /*  fixuptri.                                                                */
  4746. /*                                                                           */
  4747. /*  This routine also needs to make decisions regarding the "stacking" of    */
  4748. /*  triangles.  (Read the description of constrainededge() below before      */
  4749. /*  reading on here, so you understand the algorithm.)  If the position of   */
  4750. /*  the new point (the origin of fixuptri) indicates that the vertex before  */
  4751. /*  it on the polygon is a reflex vertex, then "stack" the triangle by       */
  4752. /*  doing nothing.  (fixuptri is an inverted triangle, which is how stacked  */
  4753. /*  triangles are identified.)                                               */
  4754. /*                                                                           */
  4755. /*  Otherwise, check whether the vertex before that was a reflex vertex.     */
  4756. /*  If so, perform an edge flip, thereby eliminating an inverted triangle    */
  4757. /*  (popping it off the stack).  The edge flip may result in the creation    */
  4758. /*  of a new inverted triangle, depending on whether or not the new vertex   */
  4759. /*  is visible to the vertex three edges behind on the polygon.              */
  4760. /*                                                                           */
  4761. /*  If neither of the two vertices behind the new vertex are reflex          */
  4762. /*  vertices, fixuptri and fartri, the triangle opposite it, are not         */
  4763. /*  inverted; hence, ensure that the edge between them is locally Delaunay.  */
  4764. /*                                                                           */
  4765. /*  `leftside' indicates whether or not fixuptri is to the left of the       */
  4766. /*  segment being inserted.  (Imagine that the segment is pointing up from   */
  4767. /*  endpoint1 to endpoint2.)                                                 */
  4768. /*                                                                           */
  4769. /*****************************************************************************/
  4770.  
  4771. void delaunayfixup(
  4772. struct triedge *fixuptri,
  4773. int leftside)
  4774. {
  4775.   struct triedge neartri;
  4776.   struct triedge fartri;
  4777.   struct edge faredge;
  4778.   point nearpoint, leftpoint, rightpoint, farpoint;
  4779.   triangle ptr;                         /* Temporary variable used by sym(). */
  4780.   shelle sptr;                      /* Temporary variable used by tspivot(). */
  4781.  
  4782.   lnext(*fixuptri, neartri);
  4783.   sym(neartri, fartri);
  4784.   /* Check if the edge opposite the origin of fixuptri can be flipped. */
  4785.   if (fartri.tri == dummytri) {
  4786.     return;
  4787.   }
  4788.   tspivot(neartri, faredge);
  4789.   if (faredge.sh != dummysh) {
  4790.     return;
  4791.   }
  4792.   /* Find all the relevant vertices. */
  4793.   apex(neartri, nearpoint);
  4794.   org(neartri, leftpoint);
  4795.   dest(neartri, rightpoint);
  4796.   apex(fartri, farpoint);
  4797.   /* Check whether the previous polygon vertex is a reflex vertex. */
  4798.   if (leftside) {
  4799.     if (counterclockwise(nearpoint, leftpoint, farpoint) <= 0.0) {
  4800.       /* leftpoint is a reflex vertex too.  Nothing can */
  4801.       /*   be done until a convex section is found.     */
  4802.       return;
  4803.     }
  4804.   } else {
  4805.     if (counterclockwise(farpoint, rightpoint, nearpoint) <= 0.0) {
  4806.       /* rightpoint is a reflex vertex too.  Nothing can */
  4807.       /*   be done until a convex section is found.      */
  4808.       return;
  4809.     }
  4810.   }
  4811.   if (counterclockwise(rightpoint, leftpoint, farpoint) > 0.0) {
  4812.     /* fartri is not an inverted triangle, and farpoint is not a reflex */
  4813.     /*   vertex.  As there are no reflex vertices, fixuptri isn't an    */
  4814.     /*   inverted triangle, either.  Hence, test the edge between the   */
  4815.     /*   triangles to ensure it is locally Delaunay.                    */
  4816.     if (incircle(leftpoint, farpoint, rightpoint, nearpoint) <= 0.0) {
  4817.       return;
  4818.     }
  4819.     /* Not locally Delaunay; go on to an edge flip. */
  4820.   }        /* else fartri is inverted; remove it from the stack by flipping. */
  4821.   flip(&neartri);
  4822.   lprevself(*fixuptri);    /* Restore the origin of fixuptri after the flip. */
  4823.   /* Recursively process the two triangles that result from the flip. */
  4824.   delaunayfixup(fixuptri, leftside);
  4825.   delaunayfixup(&fartri, leftside);
  4826. }
  4827.  
  4828. /*****************************************************************************/
  4829. /*                                                                           */
  4830. /*  constrainededge()   Force a segment into a constrained Delaunay          */
  4831. /*                      triangulation by deleting the triangles it           */
  4832. /*                      intersects, and triangulating the polygons that      */
  4833. /*                      form on each side of it.                             */
  4834. /*                                                                           */
  4835. /*  Generates a single edge connecting `endpoint1' to `endpoint2'.  The      */
  4836. /*  triangle `starttri' has `endpoint1' as its origin.  `newmark' is the     */
  4837. /*  boundary marker of the segment.                                          */
  4838. /*                                                                           */
  4839. /*  To insert a segment, every triangle whose interior intersects the        */
  4840. /*  segment is deleted.  The union of these deleted triangles is a polygon   */
  4841. /*  (which is not necessarily monotone, but is close enough), which is       */
  4842. /*  divided into two polygons by the new segment.  This routine's task is    */
  4843. /*  to generate the Delaunay triangulation of these two polygons.            */
  4844. /*                                                                           */
  4845. /*  You might think of this routine's behavior as a two-step process.  The   */
  4846. /*  first step is to walk from endpoint1 to endpoint2, flipping each edge    */
  4847. /*  encountered.  This step creates a fan of edges connected to endpoint1,   */
  4848. /*  including the desired edge to endpoint2.  The second step enforces the   */
  4849. /*  Delaunay condition on each side of the segment in an incremental manner: */
  4850. /*  proceeding along the polygon from endpoint1 to endpoint2 (this is done   */
  4851. /*  independently on each side of the segment), each vertex is "enforced"    */
  4852. /*  as if it had just been inserted, but affecting only the previous         */
  4853. /*  vertices.  The result is the same as if the vertices had been inserted   */
  4854. /*  in the order they appear on the polygon, so the result is Delaunay.      */
  4855. /*                                                                           */
  4856. /*  In truth, constrainededge() interleaves these two steps.  The procedure  */
  4857. /*  walks from endpoint1 to endpoint2, and each time an edge is encountered  */
  4858. /*  and flipped, the newly exposed vertex (at the far end of the flipped     */
  4859. /*  edge) is "enforced" upon the previously flipped edges, usually affecting */
  4860. /*  only one side of the polygon (depending upon which side of the segment   */
  4861. /*  the vertex falls on).                                                    */
  4862. /*                                                                           */
  4863. /*  The algorithm is complicated by the need to handle polygons that are not */
  4864. /*  convex.  Although the polygon is not necessarily monotone, it can be     */
  4865. /*  triangulated in a manner similar to the stack-based algorithms for       */
  4866. /*  monotone polygons.  For each reflex vertex (local concavity) of the      */
  4867. /*  polygon, there will be an inverted triangle formed by one of the edge    */
  4868. /*  flips.  (An inverted triangle is one with negative area - that is, its   */
  4869. /*  vertices are arranged in clockwise order - and is best thought of as a   */
  4870. /*  wrinkle in the fabric of the mesh.)  Each inverted triangle can be       */
  4871. /*  thought of as a reflex vertex pushed on the stack, waiting to be fixed   */
  4872. /*  later.                                                                   */
  4873. /*                                                                           */
  4874. /*  A reflex vertex is popped from the stack when a vertex is inserted that  */
  4875. /*  is visible to the reflex vertex.  (However, if the vertex behind the     */
  4876. /*  reflex vertex is not visible to the reflex vertex, a new inverted        */
  4877. /*  triangle will take its place on the stack.)  These details are handled   */
  4878. /*  by the delaunayfixup() routine above.                                    */
  4879. /*                                                                           */
  4880. /*****************************************************************************/
  4881.  
  4882. void constrainededge(
  4883. struct triedge *starttri,
  4884. point endpoint2,
  4885. int newmark)
  4886. {
  4887.   struct triedge fixuptri, fixuptri2;
  4888.   struct edge fixupedge;
  4889.   point endpoint1;
  4890.   point farpoint;
  4891.   double area;
  4892.   int collision;
  4893.   int done;
  4894.   triangle ptr;             /* Temporary variable used by sym() and oprev(). */
  4895.   shelle sptr;                      /* Temporary variable used by tspivot(). */
  4896.  
  4897.   org(*starttri, endpoint1);
  4898.   lnext(*starttri, fixuptri);
  4899.   flip(&fixuptri);
  4900.   /* `collision' indicates whether we have found a point directly */
  4901.   /*   between endpoint1 and endpoint2.                           */
  4902.   collision = 0;
  4903.   done = 0;
  4904.   do {
  4905.     org(fixuptri, farpoint);
  4906.     /* `farpoint' is the extreme point of the polygon we are "digging" */
  4907.     /*   to get from endpoint1 to endpoint2.                           */
  4908.     if ((farpoint[0] == endpoint2[0]) && (farpoint[1] == endpoint2[1])) {
  4909.       oprev(fixuptri, fixuptri2);
  4910.       /* Enforce the Delaunay condition around endpoint2. */
  4911.       delaunayfixup(&fixuptri, 0);
  4912.       delaunayfixup(&fixuptri2, 1);
  4913.       done = 1;
  4914.     } else {
  4915.       /* Check whether farpoint is to the left or right of the segment */
  4916.       /*   being inserted, to decide which edge of fixuptri to dig     */
  4917.       /*   through next.                                               */
  4918.       area = counterclockwise(endpoint1, endpoint2, farpoint);
  4919.       if (area == 0.0) {
  4920.         /* We've collided with a point between endpoint1 and endpoint2. */
  4921.         collision = 1;
  4922.         oprev(fixuptri, fixuptri2);
  4923.         /* Enforce the Delaunay condition around farpoint. */
  4924.         delaunayfixup(&fixuptri, 0);
  4925.         delaunayfixup(&fixuptri2, 1);
  4926.         done = 1;
  4927.       } else {
  4928.         if (area > 0.0) {         /* farpoint is to the left of the segment. */
  4929.           oprev(fixuptri, fixuptri2);
  4930.           /* Enforce the Delaunay condition around farpoint, on the */
  4931.           /*   left side of the segment only.                       */
  4932.           delaunayfixup(&fixuptri2, 1);
  4933.           /* Flip the edge that crosses the segment.  After the edge is */
  4934.           /*   flipped, one of its endpoints is the fan vertex, and the */
  4935.           /*   destination of fixuptri is the fan vertex.               */
  4936.           lprevself(fixuptri);
  4937.         } else {                 /* farpoint is to the right of the segment. */
  4938.           delaunayfixup(&fixuptri, 0);
  4939.           /* Flip the edge that crosses the segment.  After the edge is */
  4940.           /*   flipped, one of its endpoints is the fan vertex, and the */
  4941.           /*   destination of fixuptri is the fan vertex.               */
  4942.           oprevself(fixuptri);
  4943.         }
  4944.         /* Check for two intersecting segments. */
  4945.         tspivot(fixuptri, fixupedge);
  4946.         if (fixupedge.sh == dummysh) {
  4947.           flip(&fixuptri);   /* May create an inverted triangle on the left. */
  4948.         } else {
  4949.           /* We've collided with a segment between endpoint1 and endpoint2. */
  4950.           collision = 1;
  4951.           /* Insert a point at the intersection. */
  4952.           segmentintersection(&fixuptri, &fixupedge, endpoint2);
  4953.           done = 1;
  4954.         }
  4955.       }
  4956.     }
  4957.   } while (!done);
  4958.   /* Insert a shell edge to make the segment permanent. */
  4959.   insertshelle(&fixuptri, newmark);
  4960.   /* If there was a collision with an interceding vertex, install another */
  4961.   /*   segment connecting that vertex with endpoint2.                     */
  4962.   if (collision) {
  4963.     /* Insert the remainder of the segment. */
  4964.     if (!scoutsegment(&fixuptri, endpoint2, newmark)) {
  4965.       constrainededge(&fixuptri, endpoint2, newmark);
  4966.     }
  4967.   }
  4968. }
  4969.  
  4970. /*****************************************************************************/
  4971. /*                                                                           */
  4972. /*  insertsegment()   Insert a PSLG segment into a triangulation.            */
  4973. /*                                                                           */
  4974. /*****************************************************************************/
  4975.  
  4976. void insertsegment(
  4977. point endpoint1,
  4978. point endpoint2,
  4979. int newmark)
  4980. {
  4981.   struct triedge searchtri1, searchtri2;
  4982.   triangle encodedtri;
  4983.   point checkpoint;
  4984.   triangle ptr;                         /* Temporary variable used by sym(). */
  4985.  
  4986.   /* Find a triangle whose origin is the segment's first endpoint. */
  4987.   checkpoint = (point) NULL;
  4988.   encodedtri = point2tri(endpoint1);
  4989.   if (encodedtri != (triangle) NULL) {
  4990.     decode(encodedtri, searchtri1);
  4991.     org(searchtri1, checkpoint);
  4992.   }
  4993.   if (checkpoint != endpoint1) {
  4994.     /* Find a boundary triangle to search from. */
  4995.     searchtri1.tri = dummytri;
  4996.     searchtri1.orient = 0;
  4997.     symself(searchtri1);
  4998.     /* Search for the segment's first endpoint by point location. */
  4999.     if (locate(endpoint1, &searchtri1) != ONVERTEX) {
  5000.       vTrace(
  5001.         "Internal error in insertsegment():  Unable to locate PSLG point");
  5002.       vTrace("  (%.12g, %.12g) in triangulation.",
  5003.              endpoint1[0], endpoint1[1]);
  5004.       internalerror();
  5005.     }
  5006.   }
  5007.   /* Remember this triangle to improve subsequent point location. */
  5008.   triedgecopy(searchtri1, recenttri);
  5009.   /* Scout the beginnings of a path from the first endpoint */
  5010.   /*   toward the second.                                   */
  5011.   if (scoutsegment(&searchtri1, endpoint2, newmark)) {
  5012.     /* The segment was easily inserted. */
  5013.     return;
  5014.   }
  5015.   /* The first endpoint may have changed if a collision with an intervening */
  5016.   /*   vertex on the segment occurred.                                      */
  5017.   org(searchtri1, endpoint1);
  5018.  
  5019.   /* Find a triangle whose origin is the segment's second endpoint. */
  5020.   checkpoint = (point) NULL;
  5021.   encodedtri = point2tri(endpoint2);
  5022.   if (encodedtri != (triangle) NULL) {
  5023.     decode(encodedtri, searchtri2);
  5024.     org(searchtri2, checkpoint);
  5025.   }
  5026.   if (checkpoint != endpoint2) {
  5027.     /* Find a boundary triangle to search from. */
  5028.     searchtri2.tri = dummytri;
  5029.     searchtri2.orient = 0;
  5030.     symself(searchtri2);
  5031.     /* Search for the segment's second endpoint by point location. */
  5032.     if (locate(endpoint2, &searchtri2) != ONVERTEX) {
  5033.       vTrace(
  5034.         "Internal error in insertsegment():  Unable to locate PSLG point");
  5035.       vTrace("  (%.12g, %.12g) in triangulation.",
  5036.              endpoint2[0], endpoint2[1]);
  5037.       internalerror();
  5038.     }
  5039.   }
  5040.   /* Remember this triangle to improve subsequent point location. */
  5041.   triedgecopy(searchtri2, recenttri);
  5042.   /* Scout the beginnings of a path from the second endpoint */
  5043.   /*   toward the first.                                     */
  5044.   if (scoutsegment(&searchtri2, endpoint1, newmark)) {
  5045.     /* The segment was easily inserted. */
  5046.     return;
  5047.   }
  5048.   /* The second endpoint may have changed if a collision with an intervening */
  5049.   /*   vertex on the segment occurred.                                       */
  5050.   org(searchtri2, endpoint2);
  5051.  
  5052.     /* Insert the segment directly into the triangulation. */
  5053.     constrainededge(&searchtri1, endpoint2, newmark);
  5054. }
  5055.  
  5056. /*****************************************************************************/
  5057. /*                                                                           */
  5058. /*  formskeleton()   Create the shell edges of a triangulation, including    */
  5059. /*                   PSLG edges and edges on the convex hull.                */
  5060. /*                                                                           */
  5061. /*  The PSLG edges are read from a .poly file.  The return value is the      */
  5062. /*  number of segments in the file.                                          */
  5063. /*                                                                           */
  5064. /*****************************************************************************/
  5065.  
  5066. int formskeleton(
  5067. int *segmentlist,
  5068. int *segmentmarkerlist,
  5069. int numberofsegments)
  5070. {
  5071.   char polyfilename[6];
  5072.   int index;
  5073.   point endpoint1, endpoint2;
  5074.   int segments;
  5075.   int segmentmarkers;
  5076.   int end1, end2;
  5077.   int boundmarker;
  5078.   int i;
  5079.  
  5080.     strcpy(polyfilename, "input");
  5081.     segments = numberofsegments;
  5082.     segmentmarkers = segmentmarkerlist != (int *) NULL;
  5083.     index = 0;
  5084.     /* If segments are to be inserted, compute a mapping */
  5085.     /*   from points to triangles.                       */
  5086.     if (segments > 0) {
  5087.       makepointmap();
  5088.  
  5089.     boundmarker = 0;
  5090.     /* Read and insert the segments. */
  5091.     for (i = 1; i <= segments; i++) {
  5092.       end1 = segmentlist[index++];
  5093.       end2 = segmentlist[index++];
  5094.       if (segmentmarkers) {
  5095.         boundmarker = segmentmarkerlist[i - 1];
  5096.       }
  5097.       if ((end1 < 0) || (end1 >= inpoints)) {
  5098.       } else if ((end2 < 0) || (end2 >= inpoints)) {
  5099.       } else {
  5100.         endpoint1 = getpoint(end1);
  5101.         endpoint2 = getpoint(end2);
  5102.         if ((endpoint1[0] == endpoint2[0]) && (endpoint1[1] == endpoint2[1])) {
  5103.         } else {
  5104.           insertsegment(endpoint1, endpoint2, boundmarker);
  5105.         }
  5106.       }
  5107.     }
  5108.   } else {
  5109.     segments = 0;
  5110.   }
  5111.  
  5112.   return segments;
  5113. }
  5114.  
  5115. /**                                                                         **/
  5116. /**                                                                         **/
  5117. /********* Segment (shell edge) insertion ends here                  *********/
  5118.  
  5119. /********* Carving out holes and concavities begins here             *********/
  5120. /**                                                                         **/
  5121. /**                                                                         **/
  5122.  
  5123. /*****************************************************************************/
  5124. /*                                                                           */
  5125. /*  infecthull()   Virally infect all of the triangles of the convex hull    */
  5126. /*                 that are not protected by shell edges.  Where there are   */
  5127. /*                 shell edges, set boundary markers as appropriate.         */
  5128. /*                                                                           */
  5129. /*****************************************************************************/
  5130.  
  5131. void infecthull(void)
  5132. {
  5133.   struct triedge hulltri;
  5134.   struct triedge nexttri;
  5135.   struct triedge starttri;
  5136.   struct edge hulledge;
  5137.   triangle **deadtri;
  5138.   point horg, hdest;
  5139.   triangle ptr;                         /* Temporary variable used by sym(). */
  5140.   shelle sptr;                      /* Temporary variable used by tspivot(). */
  5141.  
  5142.   /* Find a triangle handle on the hull. */
  5143.   hulltri.tri = dummytri;
  5144.   hulltri.orient = 0;
  5145.   symself(hulltri);
  5146.   /* Remember where we started so we know when to stop. */
  5147.   triedgecopy(hulltri, starttri);
  5148.   /* Go once counterclockwise around the convex hull. */
  5149.   do {
  5150.     /* Ignore triangles that are already infected. */
  5151.     if (!infected(hulltri)) {
  5152.       /* Is the triangle protected by a shell edge? */
  5153.       tspivot(hulltri, hulledge);
  5154.       if (hulledge.sh == dummysh) {
  5155.         /* The triangle is not protected; infect it. */
  5156.         infect(hulltri);
  5157.         deadtri = (triangle **) poolalloc(&viri);
  5158.         *deadtri = hulltri.tri;
  5159.       } else {
  5160.         /* The triangle is protected; set boundary markers if appropriate. */
  5161.         if (mark(hulledge) == 0) {
  5162.           setmark(hulledge, 1);
  5163.           org(hulltri, horg);
  5164.           dest(hulltri, hdest);
  5165.           if (pointmark(horg) == 0) {
  5166.             setpointmark(horg, 1);
  5167.           }
  5168.           if (pointmark(hdest) == 0) {
  5169.             setpointmark(hdest, 1);
  5170.           }
  5171.         }
  5172.       }
  5173.     }
  5174.     /* To find the next hull edge, go clockwise around the next vertex. */
  5175.     lnextself(hulltri);
  5176.     oprev(hulltri, nexttri);
  5177.     while (nexttri.tri != dummytri) {
  5178.       triedgecopy(nexttri, hulltri);
  5179.       oprev(hulltri, nexttri);
  5180.     }
  5181.   } while (!triedgeequal(hulltri, starttri));
  5182. }
  5183.  
  5184. /*****************************************************************************/
  5185. /*                                                                           */
  5186. /*  plague()   Spread the virus from all infected triangles to any neighbors */
  5187. /*             not protected by shell edges.  Delete all infected triangles. */
  5188. /*                                                                           */
  5189. /*  This is the procedure that actually creates holes and concavities.       */
  5190. /*                                                                           */
  5191. /*  This procedure operates in two phases.  The first phase identifies all   */
  5192. /*  the triangles that will die, and marks them as infected.  They are       */
  5193. /*  marked to ensure that each triangle is added to the virus pool only      */
  5194. /*  once, so the procedure will terminate.                                   */
  5195. /*                                                                           */
  5196. /*  The second phase actually eliminates the infected triangles.  It also    */
  5197. /*  eliminates orphaned points.                                              */
  5198. /*                                                                           */
  5199. /*****************************************************************************/
  5200.  
  5201. void plague(void)
  5202. {
  5203.   struct triedge testtri;
  5204.   struct triedge neighbor;
  5205.   triangle **virusloop;
  5206.   triangle **deadtri;
  5207.   struct edge neighborshelle;
  5208.   point testpoint;
  5209.   point norg, ndest;
  5210.   int killorg;
  5211.   triangle ptr;             /* Temporary variable used by sym() and onext(). */
  5212.   shelle sptr;                      /* Temporary variable used by tspivot(). */
  5213.  
  5214.   /* Loop through all the infected triangles, spreading the virus to */
  5215.   /*   their neighbors, then to their neighbors' neighbors.          */
  5216.   traversalinit(&viri);
  5217.   virusloop = (triangle **) traverse(&viri);
  5218.   while (virusloop != (triangle **) NULL) {
  5219.     testtri.tri = *virusloop;
  5220.     /* A triangle is marked as infected by messing with one of its shell */
  5221.     /*   edges, setting it to an illegal value.  Hence, we have to       */
  5222.     /*   temporarily uninfect this triangle so that we can examine its   */
  5223.     /*   adjacent shell edges.                                           */
  5224.     uninfect(testtri);
  5225.     /* Check each of the triangle's three neighbors. */
  5226.     for (testtri.orient = 0; testtri.orient < 3; testtri.orient++) {
  5227.       /* Find the neighbor. */
  5228.       sym(testtri, neighbor);
  5229.       /* Check for a shell between the triangle and its neighbor. */
  5230.       tspivot(testtri, neighborshelle);
  5231.       /* Check if the neighbor is nonexistent or already infected. */
  5232.       if ((neighbor.tri == dummytri) || infected(neighbor)) {
  5233.         if (neighborshelle.sh != dummysh) {
  5234.           /* There is a shell edge separating the triangle from its */
  5235.           /*   neighbor, but both triangles are dying, so the shell */
  5236.           /*   edge dies too.                                       */
  5237.           shelledealloc(neighborshelle.sh);
  5238.           if (neighbor.tri != dummytri) {
  5239.             /* Make sure the shell edge doesn't get deallocated again */
  5240.             /*   later when the infected neighbor is visited.         */
  5241.             uninfect(neighbor);
  5242.             tsdissolve(neighbor);
  5243.             infect(neighbor);
  5244.           }
  5245.         }
  5246.       } else {                   /* The neighbor exists and is not infected. */
  5247.         if (neighborshelle.sh == dummysh) {
  5248.           /* There is no shell edge protecting the neighbor, so */
  5249.           /*   the neighbor becomes infected.                   */
  5250.           infect(neighbor);
  5251.           /* Ensure that the neighbor's neighbors will be infected. */
  5252.           deadtri = (triangle **) poolalloc(&viri);
  5253.           *deadtri = neighbor.tri;
  5254.         } else {               /* The neighbor is protected by a shell edge. */
  5255.           /* Remove this triangle from the shell edge. */
  5256.           stdissolve(neighborshelle);
  5257.           /* The shell edge becomes a boundary.  Set markers accordingly. */
  5258.           if (mark(neighborshelle) == 0) {
  5259.             setmark(neighborshelle, 1);
  5260.           }
  5261.           org(neighbor, norg);
  5262.           dest(neighbor, ndest);
  5263.           if (pointmark(norg) == 0) {
  5264.             setpointmark(norg, 1);
  5265.           }
  5266.           if (pointmark(ndest) == 0) {
  5267.             setpointmark(ndest, 1);
  5268.           }
  5269.         }
  5270.       }
  5271.     }
  5272.     /* Remark the triangle as infected, so it doesn't get added to the */
  5273.     /*   virus pool again.                                             */
  5274.     infect(testtri);
  5275.     virusloop = (triangle **) traverse(&viri);
  5276.   }
  5277.  
  5278.   traversalinit(&viri);
  5279.   virusloop = (triangle **) traverse(&viri);
  5280.   while (virusloop != (triangle **) NULL) {
  5281.     testtri.tri = *virusloop;
  5282.  
  5283.     /* Check each of the three corners of the triangle for elimination. */
  5284.     /*   This is done by walking around each point, checking if it is   */
  5285.     /*   still connected to at least one live triangle.                 */
  5286.     for (testtri.orient = 0; testtri.orient < 3; testtri.orient++) {
  5287.       org(testtri, testpoint);
  5288.       /* Check if the point has already been tested. */
  5289.       if (testpoint != (point) NULL) {
  5290.         killorg = 1;
  5291.         /* Mark the corner of the triangle as having been tested. */
  5292.         setorg(testtri, NULL);
  5293.         /* Walk counterclockwise about the point. */
  5294.         onext(testtri, neighbor);
  5295.         /* Stop upon reaching a boundary or the starting triangle. */
  5296.         while ((neighbor.tri != dummytri)
  5297.                && (!triedgeequal(neighbor, testtri))) {
  5298.           if (infected(neighbor)) {
  5299.             /* Mark the corner of this triangle as having been tested. */
  5300.             setorg(neighbor, NULL);
  5301.           } else {
  5302.             /* A live triangle.  The point survives. */
  5303.             killorg = 0;
  5304.           }
  5305.           /* Walk counterclockwise about the point. */
  5306.           onextself(neighbor);
  5307.         }
  5308.         /* If we reached a boundary, we must walk clockwise as well. */
  5309.         if (neighbor.tri == dummytri) {
  5310.           /* Walk clockwise about the point. */
  5311.           oprev(testtri, neighbor);
  5312.           /* Stop upon reaching a boundary. */
  5313.           while (neighbor.tri != dummytri) {
  5314.             if (infected(neighbor)) {
  5315.             /* Mark the corner of this triangle as having been tested. */
  5316.               setorg(neighbor, NULL);
  5317.             } else {
  5318.               /* A live triangle.  The point survives. */
  5319.               killorg = 0;
  5320.             }
  5321.             /* Walk clockwise about the point. */
  5322.             oprevself(neighbor);
  5323.           }
  5324.         }
  5325.         if (killorg) {
  5326.           pointdealloc(testpoint);
  5327.         }
  5328.       }
  5329.     }
  5330.  
  5331.     /* Record changes in the number of boundary edges, and disconnect */
  5332.     /*   dead triangles from their neighbors.                         */
  5333.     for (testtri.orient = 0; testtri.orient < 3; testtri.orient++) {
  5334.       sym(testtri, neighbor);
  5335.       if (neighbor.tri == dummytri) {
  5336.         /* There is no neighboring triangle on this edge, so this edge    */
  5337.         /*   is a boundary edge.  This triangle is being deleted, so this */
  5338.         /*   boundary edge is deleted.                                    */
  5339.         hullsize--;
  5340.       } else {
  5341.         /* Disconnect the triangle from its neighbor. */
  5342.         dissolve(neighbor);
  5343.         /* There is a neighboring triangle on this edge, so this edge */
  5344.         /*   becomes a boundary edge when this triangle is deleted.   */
  5345.         hullsize++;
  5346.       }
  5347.     }
  5348.     /* Return the dead triangle to the pool of triangles. */
  5349.     triangledealloc(testtri.tri);
  5350.     virusloop = (triangle **) traverse(&viri);
  5351.   }
  5352.   /* Empty the virus pool. */
  5353.   poolrestart(&viri);
  5354. }
  5355.  
  5356. /*****************************************************************************/
  5357. /*                                                                           */
  5358. /*  regionplague()   Spread regional attributes and/or area constraints      */
  5359. /*                   (from a .poly file) throughout the mesh.                */
  5360. /*                                                                           */
  5361. /*  This procedure operates in two phases.  The first phase spreads an       */
  5362. /*  attribute and/or an area constraint through a (segment-bounded) region.  */
  5363. /*  The triangles are marked to ensure that each triangle is added to the    */
  5364. /*  virus pool only once, so the procedure will terminate.                   */
  5365. /*                                                                           */
  5366. /*  The second phase uninfects all infected triangles, returning them to     */
  5367. /*  normal.                                                                  */
  5368. /*                                                                           */
  5369. /*****************************************************************************/
  5370.  
  5371. void regionplague(
  5372. double attribute,
  5373. double area)
  5374. {
  5375.   struct triedge testtri;
  5376.   struct triedge neighbor;
  5377.   triangle **virusloop;
  5378.   triangle **regiontri;
  5379.   struct edge neighborshelle;
  5380.   triangle ptr;             /* Temporary variable used by sym() and onext(). */
  5381.   shelle sptr;                      /* Temporary variable used by tspivot(). */
  5382.  
  5383.   /* Loop through all the infected triangles, spreading the attribute      */
  5384.   /*   and/or area constraint to their neighbors, then to their neighbors' */
  5385.   /*   neighbors.                                                          */
  5386.   traversalinit(&viri);
  5387.   virusloop = (triangle **) traverse(&viri);
  5388.   while (virusloop != (triangle **) NULL) {
  5389.     testtri.tri = *virusloop;
  5390.     /* A triangle is marked as infected by messing with one of its shell */
  5391.     /*   edges, setting it to an illegal value.  Hence, we have to       */
  5392.     /*   temporarily uninfect this triangle so that we can examine its   */
  5393.     /*   adjacent shell edges.                                           */
  5394.     uninfect(testtri);
  5395.     /* Check each of the triangle's three neighbors. */
  5396.     for (testtri.orient = 0; testtri.orient < 3; testtri.orient++) {
  5397.       /* Find the neighbor. */
  5398.       sym(testtri, neighbor);
  5399.       /* Check for a shell between the triangle and its neighbor. */
  5400.       tspivot(testtri, neighborshelle);
  5401.       /* Make sure the neighbor exists, is not already infected, and */
  5402.       /*   isn't protected by a shell edge.                          */
  5403.       if ((neighbor.tri != dummytri) && !infected(neighbor)
  5404.           && (neighborshelle.sh == dummysh)) {
  5405.         /* Infect the neighbor. */
  5406.         infect(neighbor);
  5407.         /* Ensure that the neighbor's neighbors will be infected. */
  5408.         regiontri = (triangle **) poolalloc(&viri);
  5409.         *regiontri = neighbor.tri;
  5410.       }
  5411.     }
  5412.     /* Remark the triangle as infected, so it doesn't get added to the */
  5413.     /*   virus pool again.                                             */
  5414.     infect(testtri);
  5415.     virusloop = (triangle **) traverse(&viri);
  5416.   }
  5417.  
  5418.   /* Uninfect all triangles. */
  5419.   traversalinit(&viri);
  5420.   virusloop = (triangle **) traverse(&viri);
  5421.   while (virusloop != (triangle **) NULL) {
  5422.     testtri.tri = *virusloop;
  5423.     uninfect(testtri);
  5424.     virusloop = (triangle **) traverse(&viri);
  5425.   }
  5426.   /* Empty the virus pool. */
  5427.   poolrestart(&viri);
  5428. }
  5429.  
  5430. /*****************************************************************************/
  5431. /*                                                                           */
  5432. /*  carveholes()   Find the holes and infect them.  Find the area            */
  5433. /*                 constraints and infect them.  Infect the convex hull.     */
  5434. /*                 Spread the infection and kill triangles.  Spread the      */
  5435. /*                 area constraints.                                         */
  5436. /*                                                                           */
  5437. /*  This routine mainly calls other routines to carry out all these          */
  5438. /*  functions.                                                               */
  5439. /*                                                                           */
  5440. /*****************************************************************************/
  5441.  
  5442. void carveholes(
  5443. double *holelist,
  5444. int holes,
  5445. double *regionlist,
  5446. int regions)
  5447. {
  5448.   struct triedge searchtri;
  5449.   struct triedge *regiontris;
  5450.   triangle **holetri;
  5451.   triangle **regiontri;
  5452.   point searchorg, searchdest;
  5453.   enum locateresult intersect;
  5454.   int i;
  5455.   triangle ptr;                         /* Temporary variable used by sym(). */
  5456.  
  5457.   if (regions > 0) {
  5458.     /* Allocate storage for the triangles in which region points fall. */
  5459.     regiontris = (struct triedge *) malloc(regions * sizeof(struct triedge));
  5460.     if (regiontris == (struct triedge *) NULL) {
  5461.       vTrace("Error:  Out of memory.");
  5462.       exit(1);
  5463.     }
  5464.   }
  5465.  
  5466.     /* Initialize a pool of viri to be used for holes, concavities, */
  5467.     /*   regional attributes, and/or regional area constraints.     */
  5468.     poolinit(&viri, sizeof(triangle *), VIRUSPERBLOCK, POINTER, 0);
  5469.  
  5470.     /* Mark as infected any unprotected triangles on the boundary. */
  5471.     /*   This is one way by which concavities are created.         */
  5472.     infecthull();
  5473.  
  5474.   if (holes > 0) {
  5475.     /* Infect each triangle in which a hole lies. */
  5476.     for (i = 0; i < 2 * holes; i += 2) {
  5477.       /* Ignore holes that aren't within the bounds of the mesh. */
  5478.       if ((holelist[i] >= xmin) && (holelist[i] <= xmax)
  5479.           && (holelist[i + 1] >= ymin) && (holelist[i + 1] <= ymax)) {
  5480.         /* Start searching from some triangle on the outer boundary. */
  5481.         searchtri.tri = dummytri;
  5482.         searchtri.orient = 0;
  5483.         symself(searchtri);
  5484.         /* Ensure that the hole is to the left of this boundary edge; */
  5485.         /*   otherwise, locate() will falsely report that the hole    */
  5486.         /*   falls within the starting triangle.                      */
  5487.         org(searchtri, searchorg);
  5488.         dest(searchtri, searchdest);
  5489.         if (counterclockwise(searchorg, searchdest, &holelist[i]) > 0.0) {
  5490.           /* Find a triangle that contains the hole. */
  5491.           intersect = locate(&holelist[i], &searchtri);
  5492.           if ((intersect != OUTSIDE) && (!infected(searchtri))) {
  5493.             /* Infect the triangle.  This is done by marking the triangle */
  5494.             /*   as infect and including the triangle in the virus pool.  */
  5495.             infect(searchtri);
  5496.             holetri = (triangle **) poolalloc(&viri);
  5497.             *holetri = searchtri.tri;
  5498.           }
  5499.         }
  5500.       }
  5501.     }
  5502.   }
  5503.  
  5504.   /* Now, we have to find all the regions BEFORE we carve the holes, because */
  5505.   /*   locate() won't work when the triangulation is no longer convex.       */
  5506.   /*   (Incidentally, this is the reason why regional attributes and area    */
  5507.   /*   constraints can't be used when refining a preexisting mesh, which     */
  5508.   /*   might not be convex; they can only be used with a freshly             */
  5509.   /*   triangulated PSLG.)                                                   */
  5510.   if (regions > 0) {
  5511.     /* Find the starting triangle for each region. */
  5512.     for (i = 0; i < regions; i++) {
  5513.       regiontris[i].tri = dummytri;
  5514.       /* Ignore region points that aren't within the bounds of the mesh. */
  5515.       if ((regionlist[4 * i] >= xmin) && (regionlist[4 * i] <= xmax) &&
  5516.           (regionlist[4 * i + 1] >= ymin) && (regionlist[4 * i + 1] <= ymax)) {
  5517.         /* Start searching from some triangle on the outer boundary. */
  5518.         searchtri.tri = dummytri;
  5519.         searchtri.orient = 0;
  5520.         symself(searchtri);
  5521.         /* Ensure that the region point is to the left of this boundary */
  5522.         /*   edge; otherwise, locate() will falsely report that the     */
  5523.         /*   region point falls within the starting triangle.           */
  5524.         org(searchtri, searchorg);
  5525.         dest(searchtri, searchdest);
  5526.         if (counterclockwise(searchorg, searchdest, ®ionlist[4 * i]) >
  5527.             0.0) {
  5528.           /* Find a triangle that contains the region point. */
  5529.           intersect = locate(®ionlist[4 * i], &searchtri);
  5530.           if ((intersect != OUTSIDE) && (!infected(searchtri))) {
  5531.             /* Record the triangle for processing after the */
  5532.             /*   holes have been carved.                    */
  5533.             triedgecopy(searchtri, regiontris[i]);
  5534.           }
  5535.         }
  5536.       }
  5537.     }
  5538.   }
  5539.  
  5540.   if (viri.items > 0) {
  5541.     /* Carve the holes and concavities. */
  5542.     plague();
  5543.   }
  5544.   /* The virus pool should be empty now. */
  5545.  
  5546.     for (i = 0; i < regions; i++) {
  5547.       if (regiontris[i].tri != dummytri) {
  5548.         /* Make sure the triangle under consideration still exists. */
  5549.         /*   It may have been eaten by the virus.                   */
  5550.         if (regiontris[i].tri[3] != (triangle) NULL) {
  5551.           /* Put one triangle in the virus pool. */
  5552.           infect(regiontris[i]);
  5553.           regiontri = (triangle **) poolalloc(&viri);
  5554.           *regiontri = regiontris[i].tri;
  5555.           /* Apply one region's attribute and/or area constraint. */
  5556.           regionplague(regionlist[4 * i + 2], regionlist[4 * i + 3]);
  5557.           /* The virus pool should be empty now. */
  5558.         }
  5559.       }
  5560.     }
  5561.  
  5562.   /* Free up memory. */
  5563.     pooldeinit(&viri);
  5564.   if (regions > 0) {
  5565.     free(regiontris);
  5566.   }
  5567. }
  5568.  
  5569. /**                                                                         **/
  5570. /**                                                                         **/
  5571. /********* Carving out holes and concavities ends here               *********/
  5572.  
  5573. /*****************************************************************************/
  5574. /*                                                                           */
  5575. /*  highorder()   Create extra nodes for quadratic subparametric elements.   */
  5576. /*                                                                           */
  5577. /*****************************************************************************/
  5578.  
  5579. void highorder(void)
  5580. {
  5581.   struct triedge triangleloop, trisym;
  5582.   struct edge checkmark;
  5583.   point newpoint;
  5584.   point torg, tdest;
  5585.   int i;
  5586.   triangle ptr;                         /* Temporary variable used by sym(). */
  5587.   shelle sptr;                      /* Temporary variable used by tspivot(). */
  5588.  
  5589.   /* The following line ensures that dead items in the pool of nodes    */
  5590.   /*   cannot be allocated for the extra nodes associated with high     */
  5591.   /*   order elements.  This ensures that the primary nodes (at the     */
  5592.   /*   corners of elements) will occur earlier in the output files, and */
  5593.   /*   have lower indices, than the extra nodes.                        */
  5594.   points.deaditemstack = (void *) NULL;
  5595.  
  5596.   traversalinit(&triangles);
  5597.   triangleloop.tri = triangletraverse();
  5598.   /* To loop over the set of edges, loop over all triangles, and look at   */
  5599.   /*   the three edges of each triangle.  If there isn't another triangle  */
  5600.   /*   adjacent to the edge, operate on the edge.  If there is another     */
  5601.   /*   adjacent triangle, operate on the edge only if the current triangle */
  5602.   /*   has a smaller pointer than its neighbor.  This way, each edge is    */
  5603.   /*   considered only once.                                               */
  5604.   while (triangleloop.tri != (triangle *) NULL) {
  5605.     for (triangleloop.orient = 0; triangleloop.orient < 3;
  5606.          triangleloop.orient++) {
  5607.       sym(triangleloop, trisym);
  5608.       if ((triangleloop.tri < trisym.tri) || (trisym.tri == dummytri)) {
  5609.         org(triangleloop, torg);
  5610.         dest(triangleloop, tdest);
  5611.         /* Create a new node in the middle of the edge.  Interpolate */
  5612.         /*   its attributes.                                         */
  5613.         newpoint = (point) poolalloc(&points);
  5614.         for (i = 0; i < 2 + nextras; i++) {
  5615.           newpoint[i] = 0.5 * (torg[i] + tdest[i]);
  5616.         }
  5617.         /* Set the new node's marker to zero or one, depending on */
  5618.         /*   whether it lies on a boundary.                       */
  5619.         setpointmark(newpoint, trisym.tri == dummytri);
  5620.         if (useshelles) {
  5621.           tspivot(triangleloop, checkmark);
  5622.           /* If this edge is a segment, transfer the marker to the new node. */
  5623.           if (checkmark.sh != dummysh) {
  5624.             setpointmark(newpoint, mark(checkmark));
  5625.           }
  5626.         }
  5627.         /* Record the new node in the (one or two) adjacent elements. */
  5628.         triangleloop.tri[highorderindex + triangleloop.orient] =
  5629.                 (triangle) newpoint;
  5630.         if (trisym.tri != dummytri) {
  5631.           trisym.tri[highorderindex + trisym.orient] = (triangle) newpoint;
  5632.         }
  5633.       }
  5634.     }
  5635.     triangleloop.tri = triangletraverse();
  5636.   }
  5637. }
  5638.  
  5639. /*****************************************************************************/
  5640. /*                                                                           */
  5641. /*  transfernodes()   Read the points from memory.                           */
  5642. /*                                                                           */
  5643. /*****************************************************************************/
  5644.  
  5645. void transfernodes(
  5646. double *pointlist,
  5647. double *pointattriblist,
  5648. int *pointmarkerlist,
  5649. int numberofpoints,
  5650. int numberofpointattribs)
  5651. {
  5652.   point pointloop;
  5653.   double x, y;
  5654.   int i, j;
  5655.   int coordindex;
  5656.   int attribindex;
  5657.  
  5658.   inpoints = numberofpoints;
  5659.   mesh_dim = 2;
  5660.   nextras = numberofpointattribs;
  5661.   if (inpoints < 3) {
  5662.     vTrace("Error:  Input must have at least three input points.");
  5663.     exit(1);
  5664.   }
  5665.  
  5666.   initializepointpool();
  5667.  
  5668.   /* Read the points. */
  5669.   coordindex = 0;
  5670.   attribindex = 0;
  5671.   for (i = 0; i < inpoints; i++) {
  5672.     pointloop = (point) poolalloc(&points);
  5673.     /* Read the point coordinates. */
  5674.     x = pointloop[0] = pointlist[coordindex++];
  5675.     y = pointloop[1] = pointlist[coordindex++];
  5676.     /* Read the point attributes. */
  5677.     for (j = 0; j < numberofpointattribs; j++) {
  5678.       pointloop[2 + j] = pointattriblist[attribindex++];
  5679.     }
  5680.     if (pointmarkerlist != (int *) NULL) {
  5681.       /* Read a point marker. */
  5682.       setpointmark(pointloop, pointmarkerlist[i]);
  5683.     } else {
  5684.       /* If no markers are specified, they default to zero. */
  5685.       setpointmark(pointloop, 0);
  5686.     }
  5687.     x = pointloop[0];
  5688.     y = pointloop[1];
  5689.     /* Determine the smallest and largest x and y coordinates. */
  5690.     if (i == 0) {
  5691.       xmin = xmax = x;
  5692.       ymin = ymax = y;
  5693.     } else {
  5694.       xmin = (x < xmin) ? x : xmin;
  5695.       xmax = (x > xmax) ? x : xmax;
  5696.       ymin = (y < ymin) ? y : ymin;
  5697.       ymax = (y > ymax) ? y : ymax;
  5698.     }
  5699.   }
  5700.  
  5701.   /* Nonexistent x value used as a flag to mark circle events in sweepline */
  5702.   /*   Delaunay algorithm.                                                 */
  5703.   xminextreme = 10 * xmin - 9 * xmax;
  5704. }
  5705.  
  5706. /*****************************************************************************/
  5707. /*                                                                           */
  5708. /*  numbernodes()   Number the points.                                       */
  5709. /*                                                                           */
  5710. /*  Each point is assigned a marker equal to its number.                     */
  5711. /*                                                                           */
  5712. /*  Used when writenodes() is not called because no .node file is written.   */
  5713. /*                                                                           */
  5714. /*****************************************************************************/
  5715.  
  5716. void numbernodes(void)
  5717. {
  5718.   point pointloop;
  5719.   int pointnumber;
  5720.  
  5721.   traversalinit(&points);
  5722.   pointloop = pointtraverse();
  5723.   pointnumber = 0;
  5724.   while (pointloop != (point) NULL) {
  5725.     setpointmark(pointloop, pointnumber);
  5726.     pointloop = pointtraverse();
  5727.     pointnumber++;
  5728.   }
  5729. }
  5730.  
  5731. /*****************************************************************************/
  5732. /*                                                                           */
  5733. /*  writeelements()   Write the triangles to an .ele file.                   */
  5734. /*                                                                           */
  5735. /*****************************************************************************/
  5736.  
  5737. void writeelements(
  5738. int **trianglelist,
  5739. double **triangleattriblist)
  5740. {
  5741.   int *tlist;
  5742.   double *talist;
  5743.   int pointindex;
  5744.   int attribindex;
  5745.   struct triedge triangleloop;
  5746.   point p1, p2, p3;
  5747.   point mid1, mid2, mid3;
  5748.   int elementnumber;
  5749.   int i;
  5750.  
  5751.   /* Allocate memory for output triangles if necessary. */
  5752.   if (*trianglelist == (int *) NULL) {
  5753.     *trianglelist = (int *) malloc(triangles.items *
  5754.                                ((order + 1) * (order + 2) / 2) * sizeof(int));
  5755.     if (*trianglelist == (int *) NULL) {
  5756.       vTrace("Error:  Out of memory.");
  5757.       exit(1);
  5758.     }
  5759.   }
  5760.   /* Allocate memory for output triangle attributes if necessary. */
  5761.   if ((eextras > 0) && (*triangleattriblist == (double *) NULL)) {
  5762.     *triangleattriblist = (double *) malloc(triangles.items * eextras *
  5763.                                           sizeof(double));
  5764.     if (*triangleattriblist == (double *) NULL) {
  5765.       vTrace("Error:  Out of memory.");
  5766.       exit(1);
  5767.     }
  5768.   }
  5769.   tlist = *trianglelist;
  5770.   talist = *triangleattriblist;
  5771.   pointindex = 0;
  5772.   attribindex = 0;
  5773.  
  5774.   traversalinit(&triangles);
  5775.   triangleloop.tri = triangletraverse();
  5776.   triangleloop.orient = 0;
  5777.   elementnumber = 0;
  5778.   while (triangleloop.tri != (triangle *) NULL) {
  5779.     org(triangleloop, p1);
  5780.     dest(triangleloop, p2);
  5781.     apex(triangleloop, p3);
  5782.     if (order == 1) {
  5783.       tlist[pointindex++] = pointmark(p1);
  5784.       tlist[pointindex++] = pointmark(p2);
  5785.       tlist[pointindex++] = pointmark(p3);
  5786.     } else {
  5787.       mid1 = (point) triangleloop.tri[highorderindex + 1];
  5788.       mid2 = (point) triangleloop.tri[highorderindex + 2];
  5789.       mid3 = (point) triangleloop.tri[highorderindex];
  5790.       tlist[pointindex++] = pointmark(p1);
  5791.       tlist[pointindex++] = pointmark(p2);
  5792.       tlist[pointindex++] = pointmark(p3);
  5793.       tlist[pointindex++] = pointmark(mid1);
  5794.       tlist[pointindex++] = pointmark(mid2);
  5795.       tlist[pointindex++] = pointmark(mid3);
  5796.     }
  5797.  
  5798.     for (i = 0; i < eextras; i++) {
  5799.       talist[attribindex++] = elemattribute(triangleloop, i);
  5800.     }
  5801.     triangleloop.tri = triangletraverse();
  5802.     elementnumber++;
  5803.   }
  5804. }
  5805.  
  5806. /*****************************************************************************/
  5807. /*                                                                           */
  5808. /*  main() or triangulate()   Gosh, do everything.                           */
  5809. /*                                                                           */
  5810. /*  The sequence is roughly as follows.  Many of these steps can be skipped, */
  5811. /*  depending on the command line switches.                                  */
  5812. /*                                                                           */
  5813. /*  - Initialize constants and parse the command line.                       */
  5814. /*  - Read the points from a file and either                                 */
  5815. /*    - triangulate them (no -r), or                                         */
  5816. /*    - read an old mesh from files and reconstruct it (-r).                 */
  5817. /*  - Insert the PSLG segments (-p), and possibly segments on the convex     */
  5818. /*      hull (-c).                                                           */
  5819. /*  - Read the holes (-p), regional attributes (-pA), and regional area      */
  5820. /*      constraints (-pa).  Carve the holes and concavities, and spread the  */
  5821. /*      regional attributes and area constraints.                            */
  5822. /*  - Enforce the constraints on minimum angle (-q) and maximum area (-a).   */
  5823. /*      Also enforce the conforming Delaunay property (-q and -a).           */
  5824. /*  - Compute the number of edges in the resulting mesh.                     */
  5825. /*  - Promote the mesh's linear triangles to higher order elements (-o).     */
  5826. /*  - Write the output files and print the statistics.                       */
  5827. /*  - Check the consistency and Delaunay property of the mesh (-C).          */
  5828. /*                                                                           */
  5829. /*****************************************************************************/
  5830.  
  5831. void triangulate(
  5832. char *triswitches,
  5833. struct triangulateio *in,
  5834. struct triangulateio *out,
  5835. struct triangulateio *vorout)
  5836. {
  5837.   double *holearray;                                        /* Array of holes. */
  5838.   double *regionarray;   /* Array of regional attributes and area constraints. */
  5839.  
  5840.   triangleinit();
  5841.   parsecommandline(1, &triswitches);
  5842.  
  5843.   transfernodes(in->pointlist, in->pointattributelist, in->pointmarkerlist,
  5844.                 in->numberofpoints, in->numberofpointattributes);
  5845.  
  5846.   hullsize = delaunay();                          /* Triangulate the points. */
  5847.  
  5848.   /* Ensure that no point can be mistaken for a triangular bounding */
  5849.   /*   box point in insertsite().                                   */
  5850.   infpoint1 = (point) NULL;
  5851.   infpoint2 = (point) NULL;
  5852.   infpoint3 = (point) NULL;
  5853.  
  5854.   if (useshelles) {
  5855.     checksegments = 1;                  /* Segments will be introduced next. */
  5856.       /* Insert PSLG segments and/or convex hull segments. */
  5857.       insegments = formskeleton(in->segmentlist, in->segmentmarkerlist,
  5858.                                 in->numberofsegments);
  5859.   }
  5860.  
  5861.     holearray = in->holelist;
  5862.     holes = in->numberofholes;
  5863.     regionarray = in->regionlist;
  5864.     regions = in->numberofregions;
  5865.       /* Carve out holes and concavities. */
  5866.       carveholes(holearray, holes, regionarray, regions);
  5867.  
  5868.   /* Compute the number of edges. */
  5869.   edges = (3l * triangles.items + hullsize) / 2l;
  5870.  
  5871.   if (order > 1) {
  5872.     highorder();             /* Promote elements to higher polynomial order. */
  5873.   }
  5874.  
  5875.   out->numberofpoints = points.items;
  5876.   out->numberofpointattributes = nextras;
  5877.   out->numberoftriangles = triangles.items;
  5878.   out->numberofcorners = (order + 1) * (order + 2) / 2;
  5879.   out->numberoftriangleattributes = eextras;
  5880.   out->numberofedges = edges;
  5881.   if (useshelles) {
  5882.     out->numberofsegments = shelles.items;
  5883.   } else {
  5884.     out->numberofsegments = hullsize;
  5885.   }
  5886.   if (vorout != (struct triangulateio *) NULL) {
  5887.     vorout->numberofpoints = triangles.items;
  5888.     vorout->numberofpointattributes = nextras;
  5889.     vorout->numberofedges = edges;
  5890.   }
  5891.     /* If not using iteration numbers, don't write a .node file if one was */
  5892.   /*   read, because the original one would be overwritten!              */
  5893.     numbernodes();                 /* We must remember to number the points. */
  5894.     writeelements(&out->trianglelist, &out->triangleattributelist);
  5895.  
  5896.     triangledeinit();
  5897. }
  5898.